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Abstract 

In this paper we develop a probabilistic micro-scale model and use 
it to study macro-scale properties of axonal transport, the processes by 
which materials are moved in the axons of neurons. By directly mod- 
eling the smallest scale interactions, we can use recent microscopic ex- 
perimental observations to infer all the parameters of the model. Then 
using techniques from queueing theory, we can predict macroscopic 
behavior in order to investigate three important biological questions: 
(1) How homogeneous are axons at stochastic equilibrium? (2) How 
quickly can axons return to stochastic equilibrium after large local 
perturbations? (3) How inhomogeneous does deposition and turnover 
make the axon? 

1 Introduction 

In all cells, one finds that proteins, membrane-bound organelles, and other 
structures (e.g. chromosomes) are transported from place to place at speeds 
much higher than diffusion. Though these transport processes are funda- 
mental to cell function, many of the underlying mechanisms, organizational 
principles, and regulatory features remain unknown. Axonal transport is 
one of the best studied systems because the transport is basically one- 
dimensional since axons are long and narrow. There are two speeds of axonal 
transport. Fast transport goes at speeds of roughly 0.2 to 0.5 meters/day 
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[20j[24j, while slow transport goes at approximately 1 millimeter /day, the 
rate of axon growth and regeneration [5] |20j . The biology and principles of 
slow transport are not yet clear [5], but the basic mechanisms of fast ax- 
onal transport were discovered in the 1980s [2j[3j[22j[34j. The model in this 
paper refers to fast axonal transport, which we will henceforth call axonal 
transport. 

The axonal transport apparatus consists of vesicles which form reversible 
chemical bonds with motor proteins that bind reversibly to microtubules 
which run parallel to the long dimension of the axon When the vesicle- 
motor protein complex is assembled on the microtubule, the complex steps 
stochastically with step size approximately 8 nanometers for kinesin and 
dynein and 10 nanometers for myosin [6j[10j[16j|31j. The vesicles enter from 
the cell body on microtubules and then detach and reattach to the transport 
mechanism at random times. 

In this paper we propose a spatial Markov-chain compartmental model 
based on these dynamics. We will assume independence of the interactions, 
and exponential wait times between events. While we address the validity of 
these assumptions in the Discussion section, we consider this a useful "first- 
order" approximation that permits study of the dynamics from both the 
perspective of individual vesicles as well as that of the full spatial system. 
Such a model unifies some earlier modeling efforts and can accommodate 
both qualitative and quantitative experimental data observed on multiple 
scales. 

In much experimental work in the 1970s and 1980s, radio-labeled amino 
acids were put into the cell bodies continuously or for a few hours. The amino 
acids were incorporated into proteins that were packaged into vesicles and 
put on the transport system so that at later times radioactivity could be 
seen moving progressively down the axons. In the continuous infusion case, 
one would see a wave of radioactivity with a sharp but slowly spreading 
wavefront propagating at constant velocity down the axon. In the case of 
infusion for a few hours one would see at long times a slowly spreading pulse 
of radioactivity that looked normally distributed. It was to understand this 
behavior that Reed and Blum constructed PDE models for axonal transport 
[3][25j[26j. These models did not have traveling wave solutions, but the data 
certainly looked like approximate traveling waves. In [27] it was shown by 
a perturbation theory argument that, in the asymptotic limit where the 
unbinding and binding rates k\ and &2 get large, the solution approaches 
a slowly spreading traveling wave or a normal pulse. Recently, in a series 
of papers, Friedman and co-workers have introduced new PDE models and 
proved these results rigorously |12| [13"] |14| |15| . 
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Probabilistic models for axonal transport were introduced and used for 
simulations already in the 1980s [30j[32j. However, rigorous work began with 
Lawler [21] in 1995 and was continued by Brooks [I] who used a continuous 
time stochastic model to show that the distribution of an individual particle 
is a spreading Gaussian at large times. Brooks also proved tail estimates for 
the central limit theorem and used them to estimate the error from normal. 

1.1 Summary of Results 

In this paper we revise the existing probabilistic models in order to ob- 
serve randomness in the system as a whole rather than exclusively from the 
particle's point of view. Our goal is to exploit experimental observations 
at multiple scales in order to make predictions about the perceived homo- 
geneity of material along the length of the axon. To this end, we propose 
in Section [2] a continuous-time Markov chain queueing model for the axonal 
transport system and estimate the order of magnitude of the various param- 
eters. The precise values of the parameters will vary for different molecular 
motors, for various types of cargo and for various animal species. 

In Section [3] we take the individual vesicle point-of-view to recover pre- 
vious results [4] [26] and show that this model does unify and extend the 
existing probabilistic and PDE models. Using standard results from re- 
newal theory we calculate the mean velocity and the near-Gaussian wave- 
front spreading of the law of the vesicle's location. Since we assume indepen- 
dence of the particle interactions, the law of an individual is equivalent to the 
distribution of an ensemble of particles released at the same time. Therefore 
the PDE governing a spatial-continuum limit of the law of a single particle 
is the same as the PDE studied in [26] (see Section l3T3j) . 

Subsequently in Section U] we adopt the full spatial system perspective to 
quantify stochasticity along the length of the axon. We begin by calculating 
in Proposition 14. II the stationary distribution of a flow-through system that 
has sustained input from the nucleus, while particles are removed upon 
reaching the distal end. The stationary distribution has a product Poisson 
structure which allows for seamless transition between spatial scales. Via 
the coefficient of variation we give a precise characterization in Section I4T21 of 
the experimenter's qualitative perception of homogeneity at the millimeter 
scale. 

The final two sections deal with model predictions. In Section [5] we 
study the non-equilibrium dynamics of recovery. Due to the product struc- 
ture of the law of the transient dynamics, behavior is determined by the 2N- 
dimensional ODE governing the means. From this we estimate the timescale 
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of return to equilibrium as a function of the lengthscale of interest. In Sec- 
tion [6] we consider the effect of deposition of material in the cell membrane. 
We find that there should be an exponential loss of material along the length 
of the cell, and conclude that a new biological mechanism must be discov- 
ered to account for the homogeneity that is a hallmark of axonal transport 
experimental observations. 

2 The model and its parameters 

Let L be the length of the axon, divided evenly into N lateral sections each 
of length 5, equal to the step size of the motor protein. Within each section, 
we disregard any further spatial geometry and take the particles to be in 
one of two states: 

• an on-transport state with mean lateral velocity v, or 

• an off-transport state with lateral velocity 0. 

We use a 2./V-dimensional continuous time Markov chain to model the 
particle dynamics: 

• Qi(t) is the number of particles in the on-transport state in section i, 

• Pi(t) is the number of particles in the off-transport state in section i. 

Each Qi and Pi has the natural numbers N = {0, 1, 2, . . .} as its state space. 
In the simplest model, we consider the following transitions and rates, 

• Lateral transport: 

(Qi, Qi+i) — >■ {Qi - 1, Qi+i + 1) at rate rQi(t), where r = v/5; 

• Switch from on-transport to off-transport: 
(Qi,Pi) -»■ [Qi ~l,P + 1) at rate fciQ»(t); 

• Switch from off-transport to on-transport: 
(Pi, Qi) -> (Pi ~l,Qi + 1) at rate k 2 Pi(t); 

• Production of new particles: 
(Qi) (Qi + 1) at rate q r; 

• Removal of particles at distal end: 
(Qn) — > (Qn - 1) at rate rQiv(i). 
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Figure 1: Double Chain 



The lateral transport rate, r, is inversely proportional to the length scale so 
that the mean number of particles per unit length is invariant with respect 
to rescaling 5. The rate of production, qor, ensures that this mean number 
per unit scales with q^. A graph of the model is depicted in Figure 1. 

In order to ensure the Markov property, we use exponential random 
variables for the waiting times between transition events. Specifically, we 
mean that after a given event we assign a new independent random variable 
to each of the 3iV + 2 possible next events, exponentially distributed with 
the appropriate rate parameter. The system of values updates according to 
the transition associated with the minimum of these waiting times. Then we 
create a new set of exponential random variables and the process proceeds 
as before. 
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This is an exactly solvable model. The advantage of computing explicit 
formulas for quantities that can be observed in experiments is that the ex- 
perimental data can then be used to determine the parameter values in the 
model. For the characterization of the approximate wavefront speed and 
spreading in Section [3] and the homogeneity calculations in Section S] we 
need order-of-magnitude estimates for the parameters. Actual parameter 
values will certainly differ depending on the particular neural tissue and 
the particular particles being transported. However, we can get order-of- 
magnitude estimates from existing data. 

First we recall that fast transport has been observed to travel at speeds of 
0.2 to 0.5 m per day. We can assume that the average velocity of particles 
while physically bound to microtubules is roughly 1 m per day, or 10~ 6 
m/s. We have already stated the assumption that the length scale of the 
individual steps satisfies 5 ~ 10 _8 m. This implies that the rate parameter 
should be r = v/5 = 10 s . 

We now turn our attention to the on-off rates rates k\ and &2 . These can 
be determined from experimentally observed run lengths on the transport 
system. Indeed, Dixit et al. [9] show that a typical run along microtubules 
for dinein and kinesin is on the order of 10 _6 m. We can compare this with 
the theoretical run length of the model to determine off-rate k\ . Within the 
model, at each step on the transport mechanism the particle has a binary 
decision to jump laterally along the transport with probability r/(k\ + r), 
or to jump off with probability k\/(r + k\). The number of jumps along the 
transport system before jumping off is therefore geometrically distributed 
on the set {0, 1, . . .} with success probability fj^r- It follows that average 
number of steps in the run is ^- and therefore the average run length is 

x 10 -8 m. Setting this equal to the average experimental run length of 
10~ 6 from [9j, we see that ^- ~ 10 2 , implying that k\ ~ Is -1 . As we will see 
in the computation of the stationary distribution in Section the ratio of 
the expected number of particles on the track to those off the track is jj^. 
Dixit [9j found that approximately 75% of the particles were motile so this 
ratio is approximately equal to 3. Since k\ ~ ls^ 1 we see that £>2 ~ 5- 

It remains to estimate q^. We will see in Proposition 14.11 that the mean 
number of particles per unit length is (1 + j^)qo = ^qo- Of course, axons 
have a large variety of diameters and larger axons will have more vesicles per 
unit length so one expects a range of values for go- However, examination 
of a large number of electron micrographs of axonal cross-sections (see for 
example [17] . Fig. 3; [18]; [23]). which are typically 100 nm thick enables 
one to estimate the number of vesicles per 100 nm segment. This number 
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is typically in the range of 10 to 100 which implies that there are 1 to 10 
vesicles per "box" in our model. Therefore qo is in the range 0.25 to 2.5, for 
various axons. 

We remark that we are ignoring some aspects of the physics and the 
biology of axonal transport. We are not including diffusion of the vesicles 
off the track. We are treating the microtubule track as though it were 
a single continuous entity from one end of the axon to the other, when 
it fact it consists of numerous, separated, microtubule fragments. And, 
we are ignoring retrograde transport and the details of the motor proteins. 
Nevertheless, this simple model will enable us to investigate the homogeneity 
questions that are the main goal of this paper. 

3 Dynamics from the Particle Perspective 

In this section we calculate properties of the stochastic dynamics by using 
standard theorems from queuing theory. In the 6 — > limit, the law of the 
location of a single particle converges to the Green's function of a linear 
partial differential equation. This enables us to obtain, as a special case, 
the asymptotic behavior of the PDE models for axonal transport in various 
asymptotic limits. 

3.1 The active transport mode 

We first consider the simple case where the particle starts at Xq = and 
stays exclusively in active transport mode. Let X t G {0, 6,26, ... , L} be the 
lateral position of a particle at time t and let rit be the number of jumps 
made by the particle as of time t. Observe, X t = 6n t . 

Proposition 3.1. Let k\ = &2 = 0, and r = v/6 > 0, then X t ~ Pois(vt). 
In particular, the mean velocity of the particle is given by jE[X t ] = v. For 
any given t > 0, in the limit as 6 — > the position of the particle satisfies 

T^ V and -^§{ Xt - vt ) t > ^^ B ^>o 
where B is a standard Brownian motion. 

Proof. Since X t = 6n t where (nt) t>0 is a Poisson process with rate r = v/6 
the set out results follow from the law of large numbers (LLN) and the 
central limit theorem (CLT) for a Poisson process. □ 
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3.2 The on/off dynamics 

We now consider a particle which undergoes transitions from on-transport 
to off-transport state and back. Denote again by X t € {0, 5, 25, ... , L} the 
lateral position of a particle at time t and let nt be the number of lateral 
transition jumps made by the particle as of time t. Observe that the particle 
will spend only a fraction of its time in active transport and hence the lateral 
speed of the particle should be slower than before. 

Proposition 3.2. Let k\,k2 > 0, and let r = v/5 > 0, then the mean 
velocity of the particle satisfies j¥,[X t ] ^ — > k ^ k2 v. 

Proof. Consider the time r for a particle to make one step of lateral trans- 
port. Before doing so a particle performs m switches from on-transport 
to off-transport and back, where m is distributed as a Geometric( fc r +r ) 

variable. The amounts of time a particle takes for each switch r- from 
on-transport to off-transport are iid variables with Exponential(/ci+r) distri- 
bution. Likewise, the amounts of time a particle takes for each switch r f 
from off-transport back to on-transport are iid variables with Exponential^) 
distribution, and are independent of the other switching times. Then 



E( 



with E[t] = (k\ + k2)/{k2r). Let n t be the number of times a particle 
makes a step of lateral transport until time t. By the Renewal Theorem 
\~E[n t \ gppj, and the result follows from X t = Sn t and 5^^ = jj^. 

□ 

Proposition 3.3. Let k\,k2 > 0, and r = v/5 > 0. In the limit as 5 — > 
the position of a particle satisfies 



X t a.s k 2 r-(Xt k 2 \ d 2k x k 2 „ 

^ i~ v -> an " vt ^—v =^ \ / — —^vBi 

t t-^oo k\ + k 2 v t k\ + k 2 > y \k\ + k 2 ) 

where B\ ~ Normal(0, 1). 

Proof. The number of lateral transport steps (nt)t>o is a renewal chain with 
E[t] = (h + k 2 )/k 2 r and Var[r] = {{k x + k 2 ) 2 + 2k x r) jr 2 k\. By the LLN 
and CLT for Renewal chains 



th a.s i , n(V± —\ jL / Var [ T ] p 

t t=^E[r]' U E[t]/^V (E[t]) 3 1 
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and now the result follows from X t = 5n t , ^fj^y = £ fcf+fc 2 = fei+fc 2 t) ' an( ^ 

2 Var[T] (5fe 2 ^ 2k 1 k 2 v 2 2kik 2 2 f x ~ n 

(EM)3 ~ fc^T^ (fci + fc 2 ) 3 ~ (fci + k 2 f V OI ~ ' 

□ 

The fact that an individual particle starting at x = will have the distri- 
bution given by Proposition 13.31 at long times was also obtained by Brooks 
[1] who proved tail estimates on the CLT and found that the remainder is 
O(^). In our approach, one can obtain the same error estimate using the 
asymptotic analysis of renewal chains, see, for example, [7]. 

3.3 Connection to Partial Differential Equations Models. 

In order to demonstrate the connection between our model and the PDEs 
seen in [27J[26j[14], we make a formal calculation regarding convergence of 
the law of the location of a single particle in the full system. 

Let (5 at = jj and let X^it) denote the lateral position of a particle in a 
system with N boxes. Define 

qN(x,t) := P{X7v(i) = x, and the particle is on track} 
PN(x,t) := F{X]y(t) = x, and the particle is off track.} 

We suppress N in the notation. By definition of the Poisson process, for 
small values of h we have 

q(x, t + h) = q(x, t)(l — {k\ + r)h) + q(x — 5, t)rh + p(x, t)k 2 h + o(h) 
p(x, t + h) = p(x, t)(l — k 2 h) + q(x, t)k\h + o{h). 

which we may reorganize as 

\{q[x, t + h) - q(x, t)) = r[q(x - S,t) - q(x, t)] - hq(x, t) + p(x, t)k 2 + 
h h 

1 o(fi) 
-(p(x, t + h) - p(x, t)) = -k 2 p(x, t) + kiq(x, t) + — — . 

Formally, we take the Taylor expansion of q in space: q[x — S, t) = q(x, t) — 
5d x q(x,t) + o(5 2 ). Taking limits as h — > and <5 — ?• yields the system of 
PDEs 

d t q(x, t) + vd x q(x, t) = ~hq(x, t) + k 2 p(x, t) (1) 
d t p(x, t) = hq(x, t) - k 2 p(x, t). (2) 
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When k\ = = k 2 , the limiting PDE is simple linear transport: {dt + 
vd x )q(x,t) = 0. The initial condition q(x, 0) = 5o(^) corresponds to the 
density of a single particle at the origin at t = 0. The time evolution via 
simple linear transport is translation of the delta function, while the time 
evolution via the equations ([I]) and ([2]) will have a spreading profile. This is 
seen in Propositions 13.11 and 13.31 since the variance of the fluctuations goes 
to in Proposition 13. II but not in Proposition 13.31 as 5 — > 0. 

In the experiments described in the introduction one sees "approximate" 
traveling waves of radioactivity in the axons in the sense that there is a slowly 
spreading wave front moving at constant velocity away from the cell body. 
The equations (3) and (4) are linear and do not have solutions that are 
bounded traveling waves. However, it was shown by a perturbation theory 
argument in [27] [26] that as e — > the solutions of 



e(d t + vd x )q £ (x,t) = -kiq £ (x,t) + k 2 p £ (x,t), q £ (0,t) = q (3) 
ed t p E (x,t) = kiq E (x,t) - k 2 p £ (x,t). (4) 

are to leading order 

q £ (x, t) = H( X ' ^ ,t), p £ {x, t) = cH{ X ^ ,t), 

where H satisfies the heat equation 

K 2 



and 



d s H(y,s) = —d yy H(y,s), H(y,0) = X(-oo,o), 



k\v 2 2kik 2 v 2 ._. 

K = 7Z i i vv ( 5 ) 



k 1 + k 2 ' (ki + k 2 ) 3 ' 

This asymptotic form is valid for small e, that is for large k\ and k 2 . How- 
ever, if we set q(x,t) = |) and p(x,t) = p £ (f , |), then q and p satisfy 
(3) and (4), so the solutions of (3) and (4) behave like approximate travel- 
ing waves for large t and large x whether or not k± and k 2 are large. The 
results suggested by these perturbation theory arguments have been proven 
rigorously by Friedman and coworkers [12] [13] [J3] [15] . 

In the case that radiolabeled particles enter the axon only for a short 
time, < t < T, the same perturbation analysis shows that for large times 
the solution of ([3]) and (|4|), where q(0, t) = for t > T is to leading order a 
spreading Gaussian with mean and variance as given in Proposition l3.3l This 



11 



corresponds to what is seen experimentally and arises asymptotically in the 
solutions of ([3]) and Q by convolving the asymptotic solution corresponding 
to the initial condition q(x,0) = 5q(x) with the indicator function X[o,T]( x )- 

4 Dynamics from the spatial system perspective 
4.1 The spatial system in equilibrium 

We are now ready to characterize the steady state dynamics induced by 
continually adding particles from the nucleus and removing them when they 
reach the distal end of the cell. 

Proposition 4.1. If the incoming rate of particles at Qq is q$r, then the 
system has the following stationary distribution 



where the {Qi} and {Pi} are mutually independent. 

Proof. Given the value of the input rate Qo = q* the generator of the process 
(Qi,Pi) in the first section is 



A q J(q,p) = [f(q + l,p) - f{q,p)}rq* + [f(q - l,p) - f(q,p)]rq 

+ [f(q - hp + 1) - f(q,p)]hq + [f(q + i,p - l) - f(q,p)]k 2 p 



If we first use f(q,p) = Qi(t) then use f(q,p) = Pi{t) and take expectations 
we get a system of ODE's governing the change in E[Qi],E[Pi] over time 



indicating that in equilibrium E[Qi] = E[Q ], E[Pj] = E[Qi]^. For pur- 
poses that will soon be clear we first assume that Qo has a constant Pois(A*) 
distribution over time rather than just being equal to A*. Let vr(g*,g,p) = 
^Xtiq*) ® ^AqG?) ® ^Xpip) where ir\ are Pois(A) distributions with rates 



stationary distribution for the process (Qi,Pi) we need to check that 





dt 

rfE[Pi](*) 
dt 



rE[Q ] - rE[Qi(*)] - fciE[Qi(t)] + k 2 E[P 1 (t)] 



A;iE[Qi(t)] -fc 2 E[Pi(t)] 




oo oo oo 



YJ2A q J(q,p)Tr(q*,q,p) = 



q t =0 q=Q p=0 
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for any choice of function / 6 T>(A qt ). 

oo oo oo \ Q* \1 \P 

q,=0 q =0p=0 H *' H ' F ' 

= e -(A»+A Q +A P ) ^ 

q*=0 q * 1 

( E E "T^f ([/(« + !.P) - + [/(? - 1,P) - /(9,p)]rg 

+ [/(g - l,p + 1) - f(q,p)}k iq + [f(q + l,p - 1) - f{q,p)]k 2 pj) 

In the inner two sums, for a fixed value of (7*, the factor multiplying f(q,p) 
for any g,p £ N 2 comes only from terms involving {q — l,q,q + 1} and 
{p - l,p,p + 1} and equals e -( A *+ A Q+ A p) times 

A Q _1 A p Ap A^ +1 Ap A^ Ap 

-rg* — T~T r( i* + 7 ^■ ~r r (g + l ) — T-y r( i 



(q-l)l p\ * q\ p\ * (9 + 1)! p\ ?! p! 

+ ( g + l)!( P -l)! fcl(g + 1) -^ lg 

W - !) ! (P+ 1)! 9! p! 

= _ rq ^ + _^_ r (g + 1) - rq + -^y-^-h(q + 1) 

<?! p! VAq q + l q + 1 Xp 



AqP + 

A Q A p / g \ 
= -7-— 7 1 -r-rg* - rg* + A*r - rg 
g! pi \A* / 

since Aq = A* and ^ = jj^. 

Summing over g* the factor multiplying f(q,p) becomes 

-(a +a p )AqA p ^ -A»A* 9 */g \ 
e pi — 2. — £_ \ 6 — - — rg* - rg* + A*r - rq) 

ql p\ ^ qj VA* / 

q* — U 

A 9 A p 

= e -(A Q +A P )_C2_^ / _ rA+ + A+r _ rq \ — 
g! p! V / 
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We've just shown that if the distribution of input particles Qo in equi- 
librium is Pois(A*) then the stationary distribution of Q\,P\ is independent 
of the distribution of Qo and is tt\ q (g> tt\ p where Xq = K[Qq], and Xp = 
E[Qo]^i/^2- It is easy to see in the above calculation that if the input Qo was 
indeed constant go then the stationary distribution of (Pi,Qi) would again 
be tt\ q <S>tt\ p with rates Xq = qo and Xp = qok\jk>2 respectively. Now, since 
the stationary distribution of Qi is the input into (Q2, P2) process, the more 
general calculation shows that in equilibrium the stationary distribution of 
(Q2, P2) is independent of the stationary distribution for Q% and is TT\ Q <Snr\ p 
with rates Xq = E[Qi] = qo and Xp = EfQij&i/A^ = go^i/^2 again. It 
follows by induction that the stationary distributions for {(Qi,Pi)} are in- 
dependent and identically distributed as i:\ Q (g) Tt\ P , Xq = qo,Xp = qo fci/fo. 
This is also an example of a clustering process satisfying the detailed balance 
conditions with linear rates discussed in Sec. 8.2 of [19] . 

□ 

4.2 Homogeneity of the axons at equilibrium 

Recall that 5 = Wnm, roughly the step size of motor proteins, and that 
axons can be up to one meter in length. Thus we are interested in phenomena 
on all the length scales 10" 5, where v = 1, 2, ... 8. Let A = 10^5; we want to 
determine how similar different segments of the axon of size A are. Let Qa 
and Pa denote the numbers of on-track and off-track particles in a segment 
of length A. 

In equilibrium, Qa and Pa are both sums of 10 v independent Poisson 
random variables with parameters Xq = qo and Xp = ^qo, respectively. 
Therefore the distributions of Qa and Pa are Poisson with parameters 
IO^Aq and lO^Ap, respectively. The mean and the variance of the num- 
ber of particles in the segment of length A is 10^(Aq + Xp). To see how 
homogeneous different slices of length A are, we consider the coefficient of 
variation, ca, which is the standard deviation divided by the mean. 

1 1 
Ca = — 1^^^=^^^^ = — 

^(Aq + A P )10" ^(1 + A:i/A; 2 )( ? oA10 8 

As indicated in Section [2J qo is in the range 0.25 to 2.5 in different axons. 
For illustrative purposes here, we will assume qo = 1. Since ki/k2 = 3, 
we see that the scale-dependent coefficient of variation ca = 1/(2 v 7 10 8 A). 
Therefore at the ten nanometer scale the coefficient of variation is simply 
1/2. At the micron scale ca = 1/20 and at the millimeter scale ca = 0.5 x 
10 _5/2 . The cutoff between "high variance" and "low variance" distributions 
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is usually considered to be when the coefficient of variation is near 1, so by 
this standard the axon is extremely homogeneous in its length at large scales. 

5 Approaching Equilibrium 

We have seen above that the axon is very homogeneous at stochastic equilib- 
rium on a space scale down to micrometers. One of the beautiful properties 
of the system of transport with reversible binding is that if it is locally out 
of equilibrium, the dynamics will automatically take it back to equilibrium. 
This is of fundamental importance to the biological function of the system 
because it means that the axon will automatically "repair" itself without 
central control of the repair process. How good is this mechanism? If a 
segment of the axon is far away from equilibrium, how long does it take 
to get back close to equilibrium? We study this question first for a single 
location and then use the estimates derived to scale the results to segments 
of any length. 

Proposition 5.1. Let (Q,P)(0) = (0,0) and let the constants a > and 
p G (0,1) satisfy the relationship a 2 p|Aoo| > 1 where Aqo is the equilibrium 
vector of (Q, P)(t). Then there exists a t* such that t > t*, 



Proof. We begin by noting that for any given f3 G (0,1), we may choose 
> such that for all t > i*, the vector of means A(i) := K[(Q, P)(t)] 
satisfies |A(i) — Aoo| < a/3|Aoo|. Then 



¥{\(Q,P)(t) - Aoo| > a|Aoo|} < F{\(Q,P)(t) - \(t)\ + \X(t) - AJ > a|Aoo|} 

<F{\(Q,P)(t)-\{t)\>a{l-P)\\ 0O \} 



Applying Chebyshev's Inequality, and observing that the variance of a Pois- 
son random variable is equal to its mean, we conclude that 



¥{\{Q,P)(t)- Aoo| > a|Aoo|} <p. 



In fact, the choice 




is sufficient, where 




F{|(Q,P)(t)-A 00 | >a|Aoo|}< 



Vsx[\(Q,P)(t)\] 

«2(1 - /^AooP 



2 



a,- 



2 (1-/3) 2 |A 



A(t)| 



■oc 



2 
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Since the initial condition for both P and Q are less than their respec- 
tive equilibrium values, each are monotonically increasing functions and the 
above reduces to 

n\(Q,P)(t) - Aoo| > a|Aoo|} < a2(1 _^ )2 | Aoo | 

for all t > t*. In order to satisfy the requirement that the right hand side 
must be less than p, we solve for f3 and find 



ay/p\\<x\ 

provided that ay/p\ Aoo| > 1. 

It remains to study the convergence of the mean and the appropriate 
choice of i*. The dynamics of the mean vector A(t) are given by the ODE 

\(t) = -Ai\(t) + q r ei (6) 

where ei is the unit vector (1,0) and 



The solution to (El) is 



k\+r -k 2 
—ki k 2 



A(t) = A oo + e- Alt (A(0)-A oo ) 
where Aqo = q^rA^ei = qo{l, This yields the estimate 

|A(t)-Aoo| < le-^Aool ^e-^IAool 
where a is the smaller of the eigenvalues of Ai, namely 

a = +k 2 +r - ^JJki + k 2 + r) 2 -4k 2 r). 

Noting that a > 0, i* may be chosen so that e~ at * < a/3, i.e. 

t* = a' 1 ln(l/(a/3)) 



□ 
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We now suppose that the whole axon is in statistical equilibrium except 
for a segment of length 5 x 10 u m in which we will assume that there are no 
particles either on or off the track. We want to know how long will it take 
for this segment to get back close to equilibrium. Of course, Proposition 
15.11 covered the case v = 0. We are interested in u = 1, . . . , 8. We imagine 
that the axon is broken up into 10 8-y segments of length 5 x W u m. In this 
rescaled system, the unbinding and binding rates per particle, k\ and k 2 
remain the same, as well as the mean on-transport velocity v. In order to 
retain this mean velocity, the rate of lateral stepping must be decreased to 
f = r x lO - " . 

The ODE for the mean vector of the rescaled system is given by 

j t \(t) = -IiA(i) + <7 Q rei (7) 

with 

* - ( kl X t ) • 

We note that the last term in ([7]) contains an r rather than an r. This is 
because the input rate is unchanged while the exit rate is diminished. 
The resulting equilibrium value is therefore rescaled as well, 

A oo = g ri r 1 e 1 = ^( fci ) A;2 ). 
Both components of this vector are of order 10 u , as expected. 

We will use the parameters discussed in Section [2] k\ = 1, &2 = \, v = 
10 _6 m/s,r = 10 2 s _1 . This implies f = 10 2_i/ s _1 . As in our analysis in 
Proposition 15-H the time to equilibrium, i*, is proportional to a -1 where a 
satisfies: 

a = ^{k 1 + k 2 + f - y/(ki +k 2 + f) 2 - 4k 2 f^ 
2k 2 f 

h + k 2 + r + v^fci + k 2 +f) 2 - Ak 2 f' 

We are most interested in the cases v > 3 (the segment has length > 10 mi- 
crons). Then k 2 r is small compared to k\. Ignoring constants and restricting 
our attention to the leading order terms gives 
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Thus, 

u ~ io"- 2 s. 

Thus, for a 10 micron segment (y = 3) the time to equilibrium is about 
10 seconds and for a 1 millimeter segment {y = 5) the time to equilibrium is 
1000 seconds or 15 minutes. Of course, the time to get close to equilibrium 
depends on the parameter p that represents what we mean by "close." 



6 Transport with Deposition 

One major feature of the biology we have ignored is that some particles in 
the off-transport state are actually deposed permanently in the membrane 
or are used for some other purpose. Such a phenomenon is easy to add to 
the model proposed above, but we must address a new qualitative feature 
of the results. Namely there is an exponential loss of material as we move 
toward the distal end. 




I 1/r 



D 2 



I 1/r 



#3 



I 1/r 



Figure 2: Fast Axonal Transport with Deposition 
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In addition to deposition, we recognize that particles are only useful in 
the membrane for a random amount of time with a given half- life, after which 
they decompose and get transported back to the cell body. Neglecting this 
retrograde transport, we augment the previous model as follows: at each 5- 
length lateral section let Di(t) be the number of particles deposed in section 
i, with N = {0, 1, 2, . . .} as the state space for Di. In addition to the previous 
dynamics we also define the following transition rates 

• (P i: Di) {Pi — 1, Dj + 1) (deposition) at rate k 3 Pi(t) 

• Di — > Di — 1 (decay) at rate 1/r where r > 

The constant r in the decay rate is the average lifetime of a deposed particle. 



6.1 Dynamics from the particle perspective 

As before, we can use the renewal theorem to determine the average speed 
of particles that are not deposed, as well as the spreading of the wavefront. 

Proposition 6.1. Let k\,k 2 ,k 3 > 0, r < oo and let r = v/S > 0, then the 
mean velocity of the particle that is not deposed into the membrane by time 
t satisfies HM^t^^v. 

Moreover, in the limit as 5 — > the position of a particle satisfies 



X t a .s k 2 + k 3 r jX t k 2 + k 3 \ ^ / 2h(k 2 + k 3 ) 

t t-s>oo k\ + k 2 + k 3 



\t k 1 + k 2 + k 3 V ) t-Kx>y (k x + k 2 + k 3 f V 1 



where B\ ~ Normal(0, 1). 

Remark 1 . Notice the mean velocity of particles is higher than when there 
is no deposition 

-v > j- — 2 v, with equality iff k 3 = 



(fci + k 2 + fc 3 ) " (h + k 2 ) 

A simple explanation for this is that particles that are still in lateral trans- 
port had returned to on-transport back from off-transport before managing 
to be deposed into the membrane, meaning that their time off-transport was 
shorter than it would have been if they did not have to beat the exponential 
clock for deposition. 



Proof. When there is deposition then by time t each particle has either 
already deposed, or is still performing lateral transport. A particle that is 
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deposed into the membrane has velocity zero from the time of its deposition 
on. A particle that can still be laterally transported must have avoided 
deposition until time t. 

If we proceed along the lines of the proof without deposition, then the 
only difference now is that the amount of time a particle takes for each switch 
from off-transport back to on-transport are i.i.d. exponential variables with 
parameter ^2+^3 rather than k 2 . Hence the amount of time r for a particle to 
make a step from one section to another has E[r] = {k\+k 2 +k 3 ) j '((k 2 +k 3 )r), 
and by the Renewal Theorem the number of steps a particle that can still 
be laterally transported satisfies|E[n t ] ^ — > jj^j^^r. 

Also Var[r] = ((A* + k 2 + fc 3 )(A; 1 + k 2 + k 3 + 2fc 1 r))/r 2 (/c 2 + k 3 ) 2 , so the 
result now follows by the Renewal LLN and CLT for n t and the fact that 

s 1 (k 2 + k 3 )r k 2 + k 3 

- 0- : — = : —V, 



E[t] k 1 + k 2 + k 3 k x + k 2 + k 3 



and 



2 Var[r] _ 5(k 2 + k 3 )v 2k 1 {k 2 + k 3 )v 2 2k x {k 2 + k 3 ) 2 
(E[r])3 h + k 2 + k 3 + (h + k 2 + fe 3 ) 3 ~ (h + k 2 + k 3 f V ° r ~ 

□ 



6.2 Dynamics from the spatial system perspective 

Once again we compute the stationary distribution of the flow-through sys- 
tem, this time with deposition. 

Proposition 6.2. Suppose the incoming rate of particles at Qo is q$r. Then 
the stationary distribution of the system is 

( k\ \ / k\k 3 T \ 

Qi ~ Pois(q p l ), Pi ~ Pois ( — - fc 9opM , A ~ 1 ^2 + k 3 q ° P J 

where p = r(k 2 + k 3 )/((r + ki)(k 2 + fc 3 ) - fci^) anrf {Qi}, {Pi}, and {Di} 
are mutually independent. 

Proof. Since (Qi,Pi) do not depend on Di the proof for the stationary dis- 
tribution of (Qi,Pi) follows along the same lines as in the case when there 
is no deposition. We first check the stationary distribution for the process 
(Qi, Pi) in the first section. For given value of Qo = q* the generator of the 
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process (Qi,P\) is given by 

A) = [/(<3i + !. p " /(Qi. p i)W* + [f(Qi - 1, Pi) - f(Qi,Pi)]rQi 
+ [/(Qi - 1, Fx + 1) - /(Qx, Pi^xQi + [/(Qj + 1, P x - 1) - /(Qi, Pi)]fc 2 Pi 

+ [/(Qi,i\-l)-/(Qi,A)]fc3A 
The system of ODE's governing the change in E[Qi],E[Pi] is now 



= E[Q ]r - E[Qi]r - E[Qi]fci + E[Pi]fc 2 



= E[Qi]fci - E[Pi]fc 2 - ElP^k 

Hence, in equilibrium 



(r + fti)(ft 2 + fa) - kik 2 

and 

E[P 1 ] = E[Q 1 h kiqo - 



fc 2 + k 3 fc 2 + fc 3 

If we let Aq =E[Qi],Ap = E[Pi], and define the distribution tt(Qo,Qi, Pi) = 
^(Qo) ® vt Aq (Qi) (8) 7r Ap (Pi), then Y^, q ,p=o A q J{q,p)ir(q*,q,p) = can 
be verified exactly as before. Since now the mean of Q\ differs from the 
mean of Qo by a multiplicative factor p < 1, and this is the mean of the 
input into Q 2 , it follows by induction that the sequence of rates for {Qi, Pi} 
Poisson random variables decays geometrically in i by a factor p l . This 
immediately implies that in equilibrium the magnitude of the variances also 
decays geometrically along the line of progression of lateral transport. 

To find the stationary distribution for Di note that if the number of 
particles at Pj = p* the generator for Di is 

A p J(d) = [f(d+l) - f(d)]k 3 p* + [f(d- 1) - f(d)](l/r)d 

It is easily checked that if P« ~ tt\ p then ir\ D is the stationary distribution 
for Di which is independent of the distribution for P« with A^ = k 3 r\p i . 

□ 

Corollary 6.3. When in equilibrium the mean and the variance of the num- 
ber of particles along the axon decay exponentially, with a loss ratio 

E[Q i+1 ] Var[Q m ] / 6h k 3 \ 

~ exp -- — — (8) 



E[Qi] Var[Qi] *\ k 2 + k 3 v 

and the same loss ratio holds for the off-transport and the deposed particles. 
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Remark 2. Notice that the rate of loss is depends on the ratio of reaction 
rates ki/k^, which we earlier estimated to be approximately 3, as well as 

the ratio k 3 /v, and that if k±, k 2 S> &3 then the loss rate is e « . 

Proof. Recalling p = r(k 2 + k 3 )J{r{k,2 + £3) + k\k 3 ) and r = v/S, we rewrite 

Observing that for any small 5, the approximation (1 + Sx)" 1 = 1 — 5x + 
(5x) 2 — . . . « e _(5:c , we have Equation flSJ). 

The same ratio holds for the other two states as well 

E[Pi\/E[P ] = Var[Pi]/Var[P ] = E[A]/E[A)] = Var[A]/Var[A)] = ^ 

□ 



6.3 Homogeneity despite deposition 

It remains to estimate the rates k 3 and r from experimental data and then 
make a formal calculation to address the volume of material loss due to 
deposition and turnover as a function of length scale. To be precise, suppose 
we want to know the percentage loss over a length A = 10" S. Then from 
the Corollary 16.31 the fraction of material loss is approximately 

E Qi -4 P V k 2 + favJ y \ k 2 + k 3 vJ v; 

So, for example, at the meter scale, v = 8, and so the loss ratio is exp(— k\k^/v{k2-\- 
fcs))- 

We now consider the parameters of a specific example, sodium pumps 
being deposed in the cell membrane. We assume that the half life of a single 
pump is approximately a week, which means r = 6 x 10 5 s. It remains to 
estimate k%. According to [36], there are roughly 1000 sodium pumps per 
micron 2 of membrane surface area. If we take the diameter of the neuron to 
be 5 microns, then the membrane surface area for one 10 nm slice is about 
0.15 fim 2 . So there are roughly 150 sodium pumps per slice. 

It is important to note that each vesicle carries roughly 100 sodium 
pumps and so the stationary distribution for the quantity D t is on average 
l/100th of the number of sodium channels which are deposed into the mem- 
brane at the location i. Taking this into account and recalling our earlier 
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parameter estimations, k% = 1, &2 = 1/3, t> = 10 -6 and go = 1 we estimate 
the average number of sodium pumps in the membrane of the first 10 nm 
slice to satisfy 150 = 100 go^i^3 T /(&2 + ^3) which implies k% = 1 x 10~ 6 . 
From the above we see that the loss ratio at the meter scale is 



which is to say that over the length of an entire meter long axon, we expect 
the concentration of sodium pumps in the membrane near the distal end to 
be roughly a 1/25 of that found near the nucleus. It is striking to note that 
over a length 0.1 m, however, the loss ratio is e - 3 = .74. This implies that 
material loss due to deposition and turnover becomes an issue exactly over 
the range of lengths of typical axons: 0.1 m to 1 m. At the lower end of 
this range, the axon is very homogeneous. At the higher end, one predicts 
significant inhomogeniety. 

7 Discussion 

In this paper, we used a spatial Markov chain model to unify several pre- 
vious approaches to modeling fast axonal transport. After inferring the 
orders of magnitude of the parameters of the model, we have used standard 
techniques from renewal theory and linear Markov chain theory to give a 
precise description of the most striking qualitative feature of axonal trans- 
port: homogeneity along the length of the axon. Furthermore we predict the 
timescale of recovery to this homogeneous state after removing the vesicles 
from a segment as a function of the length scale of the segment. 

In creating the model, we have become aware of two features of the 
biology which have not been adequately studied: 1) Deposition and decay 
of proteins; and 2) characterization of the process by which unbound vesicles 
reattach to the transport mechanism. 

Regarding deposition and decay, we mentioned in Section [6] that the 
model unequivocally predicts a loss of material down the length of the axon. 
Furthermore, it seems that certain biological features that we have chosen 
to ignore - such as retrograde transport and a location-dependent unbinding 
rate [9] - would only serve to increase the non-homogeneity. In the case of 
sodium pumps, it is not immediately clear if there is a contradiction between 
deposition loss and homogeneity, but we do see a sensitive dependence on 
the length of the axon. Furthermore, it stands to reason that for other types 
of cargo the loss may significant. We believe that this raises an interesting 
biological issue, and mention a few possible mechanisms for overcoming the 




(10) 
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material loss: first, it is possible that proteins can be synthesized in axons 
|33j ; second, there may be some signal that prevents material from deposing 
in a region that is already "full;" and third, it may be that for some materials 
"deposition" is a reversible process. 

Regarding reattachment, we note that we have modeled all event wait 
times as exponential random variables, but this is certainly a simplification. 
As an example, the stepping process of kinesin is a well-studied though still 
not completely understood phenomenon. Much work has focused on as- 
sessing the dependence of the mean rate of translocation on both the load 
and the local concentration of ATP |35] [29J. Implicit in this analysis is 
the assumption of exponential wait times with state dependent rate param- 
eters. However, when fitting to data and matching dispersion information 
the authors in found it necessary to generalize the wait time distribu- 
tion. This was followed by more detailed models for which it was shown that 
load carrying could in fact regularize the stepping times of kinesin motors 

m m- 

Generalizing the wait times does not significantly affect the particle per- 
spective results of Sections 13.21 and 16.11 This is because renewal theory, 
which depends on independence of the wait times, is robust with respect to 
such changes. However, the spatial system results may be affected in that 
the full process is no longer Markov. While one may still expect a some- 
thing like a stationary distribution for the flow-through system (where new 
particles enter from the nucleus and particles are removed from the system 
at the distal end), the analysis of the approach to equilibrium may change. 
In the absence of direct observation of the phenomenon, we refrain from this 
more detailed analysis. 

In light of this known need for generalized wait times in the stepping 
process, it seems likely that detailed observation of the rebinding process 
will call for new models as well. Recall that when a vesicle unbinds from 
a microtubule it is unclear whether it typically rebinds to the same micro- 
tubule or if it explores the region significantly via diffusion before finding 
a different microtubule to bind to. In the latter case, a more appropriate 
model for rebinding time would be to solve some kind of first passage time 
problem and use that distribution for the rebinding wait. 
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A STOCHASTIC COMPARTMENTAL MODEL FOR FAST AXONAL 

TRANSPORT 



LEA POPOVIC*, SCOTT A. MCKINLEY" 1 " , AND MICHAEL C. REED* 

Abstract. In this paper we develop a probabilistic micro-scale compartmental model and use 
it to study macro-scale properties of axonal transport, the process by which intracellular cargo 
is moved in the axons of neurons. By directly modeling the smallest scale interactions, we can 
use recent microscopic experimental observations to infer all the parameters of the model. Then, 
using techniques from probability theory, we compute asymptotic limits of the stochastic behavior 
of individual motor-cargo complexes, while also characterizing both equilibrium and non-equilibrium 
ensemble behavior. We use these results in order to investigate three important biological questions: 
(1) How homogeneous are axons at stochastic equilibrium? (2) How quickly can axons return to 
stochastic equilibrium after large local perturbations? (3) How is our understanding of delivery time 
to a depleted target region changed by taking the whole cell point-of-view? 

1. Introduction. In all cells, one finds that proteins, membrane-bound or- 
ganelles, and other structures (e.g. chromosomes) are transported from place to place 
at speeds much higher than diffusion. Though these transport processes are funda- 
mental to cell function, many of the underlying mechanisms, organizational principles, 
and regulatory features remain unknown. Axonal transport is one of the best studied 
systems because the transport is basically one-dimensional since axons are long and 
narrow. There are two speeds of axonal transport. Fast transport goes at speeds of 
roughly 0.2 to 0.5 meters/day [27][33], while slow transport goes at approximately 1 
millimeter/day, the rate of axon growth and regeneration [B]|27|. The biology and 
principles of slow transport are not yet clear [5J, but the basic mechanisms of fast 
axonal transport were discovered in the 1980s [2] [3] [12] [13] • The model in this paper 
refers to fast axonal transport, which we will henceforth call axonal transport. 

The axonal transport apparatus consists of vesicles which form reversible chemical 
bonds with motor proteins that bind reversibly to microtubules which run parallel to 
the long dimension of the axon pQ . When the vesicle-motor protein complex is assem- 
bled on the microtubule, the complex steps stochastically with step size approximately 
8 nanometers for kinesin and dynein and 10 nanometers for myosin [7] [12 [18 40 . The 
vesicles enter from the cell body on microtubules and then detach and reattach to the 
transport mechanism at random times. 

In this paper we propose a spatial Markov-chain compartmental model based on 
these dynamics. We will assume independence of the interactions, and exponential 
wait times between events. While we address the validity of these assumptions in the 
Discussion section, we consider this a useful "first-order" approximation that permits 
study of the dynamics from both the perspective of individual vesicles as well as that 
of the full spatial system. Such a model unifies all earlier deterministic and stochastic 
modeling efforts and can accommodate both qualitative and quantitative experimental 
data observed on multiple scales. 

In much experimental work in the 1970s and 1980s, radio-labeled amino acids 
were put into the cell bodies continuously or for a few hours. The amino acids were 
incorporated into proteins that were packaged into vesicles and put on the transport 
system so that at later times radioactivity could be seen moving progressively down 
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the axons. In the continuous infusion case, one would see a wave of radioactivity 
with a sharp but slowly spreading wavefront propagating at constant velocity down 
the axon. In the case of infusion for a few hours one would see at long times a 
slowly spreading pulse of radioactivity that looked normally distributed. It was to 
understand this behavior that Reed and Blum constructed PDE models for axonal 
transport [3] [34] [35] . These models did not have traveling wave solutions, but the 
data certainly looked like approximate traveling waves. In [36] it was shown by a 
perturbation theory argument that, in the asymptotic limit where the unbinding and 
binding rates ki and k\ get large, the solution approaches a slowly spreading traveling 
wave or a normal pulse. Recently, in a series of papers, Friedman and co-workers have 
introduced new PDE models and proved these results rigorously [T4] [T5] [HI] [T7] . 

Probabilistic models for axonal transport were introduced and used for simula- 
tions already in the 1980s [39 41.. However, rigorous work began with Lawler [25] in 
1995 and was continued by Brooks [5] who used a continuous time stochastic model 
to show that the distribution of an individual particle is a spreading Gaussian at large 
times. Brooks also proved tail estimates for the central limit theorem and used them 
to estimate the error from normal. Independently, Bressloff [4] developed a discrete 
stepping model and performed an analysis under the assumption that the rate of un- 
binding and binding to transport is fast relative to lateral velocity over the length 
scale of interest. The author derived a characterization of the spreading wavefront of 
a particle entering at the nucleus and traveling to the distal end. This model served as 
the basis for later investigations by Newby and Bressloff [31] [32] wherein the authors 
characterize the axonal transport system as an intermittent search for hidden targets. 

In this paper we revise the existing probabilistic models in order to study ran- 
domness in the system as a whole rather than exclusively from the point of view of 
an individual particle. Our goal is not only to recover and generalize previous results, 
but also to investigate three specific, biologically important, properties of the whole 
stochastic system. 

1.1. Summary of Results. In Section [2] we create a continuous-time Markov 
chain queueing model for the axonal transport system. We show how to use experi- 
mental data to determine (or estimate) all the parameters of the model. 

In Section [3] we take the individual vesicle point-of-view. We prove the asymp- 
totic forms in [36] with rigorous error estimates. We show that in the limit as the 
compartment size becomes small our model becomes the probabilistic model of [5]- 
We also show that in the limit as the length of the axon and time become large (with 
the scale of axon length on the order of the squared scale of time) our model becomes 
the PDE model of [35] ■ Since we assume that particles are independent, the time evo- 
lution of the law of an individual will reflect the behavior of an ensemble of particles 
released at the same time. 

In Section [4] we adopt the full spatial system perspective to quantify stochasticity 
along the length of the axon. We begin by calculating in Proposition ^, li the stationary 
distribution of a flow-through system that has sustained input from the nucleus, while 
particles are removed upon reaching the distal end. The stationary distribution has a 
product Poisson structure which allows for seamless transition between spatial scales. 

With this mathematical model, we are able to make precise statements about 
three biologically important system properties. In the stationary distribution, the 
number of vesicles in each slice of the axon is independent and identically distributed 
(Proposition 14. ip , however this in and of itself is not sufficient to account for the 
sense that samples taken for different parts of the cell "look the same." We compute in 
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Section B~2"l the coefficient of variation for the number of vesicles in sections of different 
length and show that the coefficient of variation is low for all but the smallest length 
scales. In Section 14.31 we study the intermittent search problem posed by Newby 
and Bressloff [31 from the system point of view. Efficient transport to locations that 
need material must balance the speed of transport of material from the nucleus to the 
distal end of the cell with the rates of dissociation from the transport apparatus. We 
calculate the expected hitting time for a hidden target by all vesicles in the system. 
In so doing we encounter the counterintuitive result that while increasing the velocity 
of the motors while on transport increases the chance of any particular vesicle missing 
the target, the expected hitting time by the system actually decreases. 

This hitting time approach is natural for needed material that is sparsely dis- 
tributed throughout the axon, but when the needed cargo in question is more com- 
mon, the time to replenishment is better addressed through the ODE approach that 
we develop in Section l4~4l Due to the product structure of the law of the transient 
dynamics, this non-equilibrium behavior is determined by the 27V-dimensional ODE 
governing the means. From this we estimate the timescale of return to equilibrium as 
a function of the length scale of interest. 

2. The model and its parameters. Let L be the length of the axon, divided 
evenly into N = L/S lateral sections each of length 6, equal to the step size of the 
motor protein. Within each section, we disregard any further spatial geometry and 
take the particles to be in one of two states: 

■ an on-transport state that steps laterally at a rate r = v/S per section, or 

• an off-transport state that does not step laterally. 

We use a 2A r -dimensional continuous time Markov chain {[Qi[t), Pj(i)), i = 1, . . . , N} 
to model the particle dynamics, where 

• Qi(t) is the number of particles at time t in on-transport state in section i, 

■ Pi(t) is the number of particles at time t in off-transport state in section i. 
Definition 2.1. (Stochastic compartmental model) Let [{[Qi[t), Pi(t)),i — 1, . . . , iV})t>o 

be a continuous time Markov chain on the state space N N x N N with the following tran- 
sitions and time dependent rates: 

■ (Lateral transport) 

(Qi,Qi+i) -> [Qi - 1, Qi+i + 1) at rate rQi(t); 

■ (Switch from on-transport to off-transport) 
[Qi, Pi) -> [Qi -l,Pi + 1) at rate k 2 Qi[t); 

■ (Switch from off-transport to on-transport) 
[Pi, Qi) -> [Pi -l,Qi + 1) at rate k x Pi[t); 

■ (Production of new particles) 
(Ql) [Ql + 1) at rate rq ; 

■ (Removal of particles at distal end) 
[Qn) -> [Qn — 1) at rate rQ N (t). 

The lateral transport rate, r = v/6, is inversely proportional to the length scale 
so that the mean number of particles per unit length is invariant with respect to 
rescaling 5. We will assume that the rate of production qo — 5po, for some constant 
Po > 0, in order for the mean number of particles in each compartment to scale with 
the size 5 of each compartment. This will imply that the mean number of particles 
per unit length scales as po- A graph of the model is depicted in Figure 1. 

In order to insure the Markov property, we use exponential random variables for 
the waiting times between transition events. Specifically, we mean that after a given 
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Fig. 2.1. On- and- off Transport Chain 

event we assign a new independent random variable to each of the 3-/V+2 possible next 
events, exponentially distributed with the appropriate rate parameter. The system 
of values updates according to the transition associated with the minimum of these 
waiting times. Then we create a new set of exponential random variables and the 
process proceeds as before. 

The advantage of computing explicit formulas for quantities that can be observed 
in experiments is that the experimental data can then be used to determine the 
parameter values in the model. For the characterization of the approximate wavefront 
speed and spreading in Section [3] and the homogeneity calculations in Section [4] we 
need order-of-magnitude estimates for the parameters. Actual parameter values will 
certainly differ depending on the particular neural tissue and the particular particles 
being transported. However, we can get order-of-magnitude estimates from existing 
data. 

First we recall that fast transport has been observed to travel at speeds of 0.2 to 
0.5 m per day. We can assume that the average velocity of particles while physically 
bound to microtubules is roughly 1 m per day, or v = 10 -6 m/s. We assume that the 
compartment size scales as the length scale of the individual steps of the motor protein, 
so 5 <~*j 10~ 8 m. This implies that the rate parameter should be r = v/S = 100s _1 . 

We now turn our attention to the on-off rates rates ki and k\ . These can be de- 
termined from experimentally observed run lengths on the transport system. Indeed, 
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Dixit et al. |10) show that a typical run along microtubules for dinein and kinesin 
is on the order of 10 _6 m. We can compare this with the theoretical run length of 
the model to determine off-rate k^. Within the model, at each step on the transport 
mechanism the particle has a binary decision to jump laterally along the transport 
with probability rj (ki + r), or to jump off with probability kij[r + fe). The number 
of jumps along the transport system before jumping off is therefore geometrically dis- 
tributed on the set {0, 1, . . .} with success probability r ,^ r fc; . It follows that average 
number of steps in the run is and therefore the average run length is x 10~ 8 m. 
Setting this equal to the average experimental run length of 10 -6 from [TO], we see 
that ~ 100, implying that k^ ~ Is -1 . As we will see in the computation of the 
stationary distribution in Section 14.11 the ratio of the expected number of particles 
on the track to those off the track is Dixit [TO] found that approximately 75% of 
the particles were motile so this ratio is approximately equal to 3. Since k-i ~ Is -1 
we see that ki *~ |. 

It remains to estimate qo- We will see in Proposition 14. II that the mean number 
of particles per compartment is (1 + jr)qo = 4go- Of course, axons have a large 
variety of diameters and larger axons will have more vesicles per unit length so one 
expects a range of values for qg. However, examination of a large number of electron 
micrographs of axonal cross-sections (see for example [20], Fig. 3; [21] ; [30] ), which 
are typically 100 nm thick enables one to estimate the number of vesicles per 100 nm 
segment. This number is typically in the range of 10 to 100 which implies that there 
are 1 to 10 vesicles per compartment in our model. Therefore go is in the range 0.25 
to 2.5, for various axons. 

We remark that we are ignoring some aspects of the physics and the biology of 
axonal transport. We are not including diffusion of the vesicles off the track. We 
are treating the microtubule track as though it were a single continuous entity from 
one end of the axon to the other, when it fact it consists of numerous, separated, 
microtubule fragments. And, we are ignoring retrograde transport and the details of 
the motor proteins. Nevertheless, this simple model will enable us to investigate the 
homogeneity questions that are the main goal of this paper. 

3. Dynamics from the Particle Perspective. In this section we calculate 
properties of the stochastic dynamics by using stochastic convergence theorems and 
stochastic averaging theorems from probability theory. We first see that, in the S — > 
limit, the law of the location of a single particle corresponds to that of a particle with 
a piecewise linear Markov motion. We then show this law can be approximated by the 
Green's function of a linear partial differential equation. This enables us to obtain, 
as a special case, the asymptotic behavior of the PDE models for axonal transport in 
an either rate limiting or perturbed setting. 

3.1. The active transport mode. We first consider the simple case where 
the particle starts at X^ — and stays exclusively in active transport mode. Let 

€ {0, S, 2(5, ... , L} be the lateral position of a particle at time t and let nf be the 
number of jumps made by the particle as of time t. Observe, X s = Sn s . 

Lemma 3.1. Let ki = k\ = 0, and r = v/S > 0, then the position of the particle 
satisfies, for any time t < oo 




s<t 



a.s. 
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and, for B a standard Brownian motion 

in distribution on the Skorokhod space of cadlag (right continuous left limited) func- 
tions. 

Proof. Since n s is a Poisson process with rate r = v/S, defining N t := n| t we get 
a Poisson process N with rate v, and we have X s = SN./§. Our results then follow 
directly from the functional law of large numbers (FLLN) and the functional central 
limit theorem (FCLT) for the Poisson process N. □ 

3.2. The on/off dynamics. We now consider a particle which undergoes tran- 
sitions from on-transport to off-transport state and back. Denote again by Xf £ 
{0, S, 25, ... , L} the lateral position of a particle at time t and let nf be the number of 
lateral transition jumps made by the particle as of time t. Observe that the particle 
will spend only a fraction of its time in active transport and hence the lateral speed 
of the particle should be slower than before. 

A non-compartmental stochastic model for axonal transport, introduced by Brooks 
in (5], is as follows. A particle can be in one of two states: 

• an on-transport state with deterministic lateral velocity v, or 

• an off-transport state with lateral velocity 0. 

We use a 2-dimensional Markov process to model the particle dynamics, where 

• X t be the lateral position of this particle at time t, 

• £t be the indicator for whether it is on (1) or off (0) transport at time t. 
Definition 3.2. (Stochastic non-compartmental model) Let (Xt,£,t)t>o be a 

piecewise- linear Markov process with values in (R+, {0, 1}) started at (A"o, £o) = (0, 1) 
with the following dynamics: 

■ (Switch from on-transport to off-transport) 
(X t , I) ->• (X t ,0) at rate k x 

■ (Switch from off -transport to on-transport) 
{X t ,0)-> (X u l) ratek 2 

■ (Lateral travel) X t — J Q v£ s ds 

The path of (A 4 ) f >o consists of alternating sequence of Exponential(fc 2 ) stretches of 
time where the lateral position increases linearly with speed v, and Exponential(fci) 
stretches of time where it remains constant. 

Proposition 3.3. Let fc 2 ,fci > 0, and r = v/S > 0, then the position of the 
particle converges 

(Xt)t>o =>(X t )t>o 

(5— >0 

in distribution on the Skorokhod space of cadlag functions. 

Proof. If we let £ s be the indicator of whether the particle in the compartmental 
model is on (£* = 1) or off (£ 5 = 0) transport, then (X$,$)t>o is a strong Markov 
process. We will see that £ s is a continuous-time Markov chain on {0, 1} (independent 
of 5), and conditionally on £ s the transition law of X is easily expressed. Likewise, 
for the non-compartmental model above, (X t ,£,t)t>o is a strong Markov process, with 
£ the same continuous-time Markov chain on {0, 1} as £ s , and conditionally on £ the 
change in X is easily given in terms of its linear speed and £. 

We will start by showing that for any t > 0, (Xf, converges to (Xt,£t) in 
distribution as S — ¥ 0. We then show that the finite dimensional distributions of 
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(X s , £ <5 ) converge to those of (X, £). A tightness argument finally implies (Lemma 
16.2 and Theorem 16.3 [22]) that (X" 5 ,^) converges to (X, £) in distribution on the 
Skorokhod space of cadlag processes. 

Suppose that initially the particle is on transport at xq, so (Xq,£q) = (xq,1). 
The first time T\ = inf{i > : £f = 0} at which the particle steps off transport 
has Exponential^) distribution, irrespective of S. The first subsequent increment in 
time o i = inf{f > : £, s 1+t = 0} after which the particle steps back on transport has 
Exponential(fci) distribution, irrespective of 5 as well. This is repeated, and £ 5 is a 
simple continuous-time Markov chain on {0, 1} with transition rates £; 2 and k\, from 
1 — > and — > 1, respectively. 

Until time t\ the particle behaves as if it were in active transport (fc 2 = k\ = 0) 
and conditionally on the value of T\ , for any < s < n the lateral change in position 
over time s, Xf — Xq, is (5nf where nf is a Poisson(rs) random variable. Moreover, 
by results of Proposition 13.11 conditionally on the value of n , we have 

sup \(X%-X$)-vs\ — > a.s. 
0<s<n 

To get the unconditioned law of X* — Xq, we observe that the number of boxes 
the particle traverses before it steps off has a Geometric (fc 2 /(fc 2 + r)) distribution 
(note that a Poisson rate r process sampled at an Exponential (fc 2 ) time independent 
of the process has this distribution). Since fc 2 /<5(fc 2 + r) — » fc 2 /w, as <5 — >■ 0, X* — Xq = 
5n s Ti converges in distribution to Exponential(fc 2 /w) variable. 

Consider now the particle in the non-compartmental model. It is immediate from 
the definition of the model that £ has the same law of a continuous-time Markov chain 
on {0, 1} with transition rates fc 2 and k\, from 1 — > and — > 1, as £ A . Since we prove 
our convergence in law results by first conditioning on the values of £ s and £, we can 
without loss of generality henceforth assume £ 5 = £ a.s., and drop its superscript. 

Suppose that initially the particle in the non-compartmental model is on transport 
at xq, so (Xo, £o) = ( x 0i !)• Conditionally on t\, for all < s < Tx, X s — Xq — vs and 
X Tl — Xq = vt\, hence unconditionally, X Tl — Xo is an Exponential(fc 2 /v) random 
variable. We now have both sup 0<s<ri |(Xf — Xq) — (X s — Xq) — > a.s. and 
X s 1 — Xq X Tl — Xq. 

Between times t\ and <j\ the particle in both models stays in place, so condition- 
ally on values of ti,cti, and ofX£,X Tl , sup Ti < s < Ti+(Ti | (Xf -X*J - (X s -X Tl )\ = 0, 
and X^ i+(Ti — Xq => X n + (7l — Xo. At time t\+&\, the same process starts over from ini- 
tial values (X* , 1) and (X Tl , 1) in the compartmental and non-compartmental model, 
respectively. 

Let Ti , o\ , r 2 , . . . , be the sequence of time increments between the consecutive 
times when the particle in both models gets off and gets back on transport, let Co = 
and for i > 1 

n = inf{< > : £a z l +t = 0}, a, = inf{< > : £ Tl+4 = 1} 

Then (rj)j>i and (<r^)j>i are independent sequences of i.i.d. Exponential (fc 2 ) and 
Exponential(fci) variables, respectively. For any t > 0, let rjt be the number of times 
the particle in either model gets back on transport until time t, and rj' t the number of 
times it gets off, 



k k' fc'-l 

T] t = inf{fc > : J2(n + o-i) < *}, rf t = inf{fc' > : ]T n + ^ <n < t} 

i—1 i—1 i—1 
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Note that, rj' t = r\ t iff £ t = 1, and rj' t = r\t + 1 iff £t = 0. Let f t be the last time 
before time t that the particle changed whether it was on or off transport, that is, 
Tt — sup{0 < s < t : £ a - ^ £, s }- Then, we have 

m 

E T * + E ^ if v't = vt + !• 

i=l i=l 

If ?7j = ??t, then from time Tt to i the particle is in active transport, and the same 
convergence argument as before implies that conditionally on the values of r* and X s , 

sup \(X S S - X s Tt ) - v(s - Tt )\ -> a.s. 

T t <S<t 

Also X* — Xq = <5n£ where n£ is the number of boxes the particle traverses by time 
t. Conditionally on the value of m, n s is a sum of rj t i.i.d. Gcometric(fc2/(fc2 + f)) 
random variables, hence X° — Xq converges in distribution to a sum of r\t i.i.d. 
Exponential(fc2/w) random variables. In the non-compartmental model, if r/' t — rjt, 
then conditionally on the values of rjt and Tt, for Tt < s < t, X s — X Tt = v(s — Tt) and 
conditionally only on the value of rjt, X s — Xo is a sum of r\t i.i.d. Exponential^/ - ^) 
variables. Hence, conditionally on rj t and T f , sup rt<s<t | (X* — X s t ) — (X s — X Tt )\ — > 
a.s., and conditionally only on jft, Xl — X^ X Tt — Xo. 

If r)' t = rjt + 1, then from time T t to £ the particle is in both models stays in place, 
so conditionally on r] t and Tt, sup Tf<s<t | (X s — X S J — (X s — X Tt )| = a.s.. Also, con- 
ditionally only on the value of r/' t , n s t is a sum of r]' t i.i.d. Geometric^/ (&2 + r)) vari- 
ables, hence X s — Xq converges in distribution to a sum of rj' t i.i.d. Exponential (fe/f) 
variables. In the non-compartmental model, if r/' t = r\t + 1, conditionally only on the 
value of r]'t, X s is a sum of r\ t i.i.d. Exponential^/^) variables. Hence, conditionally 
only on r,' t , x( - X s => X T% - X . 

Now, integrating over the possible values of ryt,^ and Tt, we get that for any t > 0, 
(X s , £ s ) =>■ (X t ,£t). Convergence of finite dimensional distributions follows from an 
iterative use of the Markov property of (X 5 ,^ 5 ) and (X, £), and the fact that the 
increments of both (X s , £ s ) and (X, £) are stationary. 

In order to verify tightness, Theorem 16.11 [22], of the sequence of Markov pro- 
cesses {(X s , £ s )}s>o, because (X s , £*) has stationary increments and is strong Markov, 
it will suffice to check that for any e > 

lim limsupP{||(X^,efj - (X S ,H S o)\\ > e} = 

h^Q x^ 

where ||(a;i, — (x2, £,2)\\ — \x\ — 1 + ICi — £2! is a distance metric on R + x {0, 1}. For 
any 5 > 0, the first change in the continuous-time Markov chain happens after an 
Exponential^) time (where k = &2 or k = k\ depending on whether £q = 1 or £q = 0), 
and is independent of 5. Hence, at time h later, P{£f 7^ £o} < 1 — e~ kh . Moreover, 
irrespective of the value of £ s , at time h later the value of |X^ — Xq| < 8n s h where n s is 
a Poisson process with rate r = v/S. Hence, P{\X S - X s \ > e} < SE[n s h ]/e = u/i/e- 
Combining the two gives 

~ (^o^o)ll > e} < 1 - e- feh + wA/e for any 5 > 0, 
and the desired limit follows. □ 
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The process {X s , £_ s ) : t £ [0, oo) ^ {X s , £ s ) G SZ + x {0,1} is a Markov process 
with cadlag paths whose generator is given by 

A 5 f(x,Z)=rZ[f(x + d,0-f(x,0] 

+ k 2 t [fix, C - 1) - f{x, 0] + fci(l - [/(a:, f + 1) - /(*, 0] 

for all / e V{A S ) = C°(6Z+ x {0, 1}). 

The piecewise linear process (X, £) : i € [0, oo) H> (X t ,^) e K + x {0,1} is a 
Markov process with continuous paths whose generator is the closure of the operator 

Af(x,0 = vtd x f(x,t) 

+ k 2 £ [f{x, £ - 1) - fix, 0} + fci(l - [/(*, £ + 1) - /(x, 0] 

for all / <= 2?(yl 5 ) = C 1 '°(R + x {0, 1}). 

Letting l s : SZ + x {0, 1} H> E + x {0, 1} be an embedding, and f s — f o i s ' , then 
A s f s -> A/ as (5 -> for all / € C 1 '°(K + x {0, 1}) imply that the finite dimensional 
distributions of {X s , £ 5 ) converge to those of {X, £). Verification of additional con- 
ditions, see Theorem 19.25 of [55], would also imply convergence of processes with 
generators {{A s ,V(A s ))}s >0 to the process with generator (A, V{A 5 )) , however, we 
thought this way of showing convergence in law was not as instructive. 

The fact that an individual particle will have the distribution given by Propo- 
sition 13.31 as the size of the boxes decreases means that our model is a microscopic 
version of the stochastic model used by Brooks [5], and that the hydrodynamic limit 
of our model as S — >• is equal to the macroscopic stochastic model from [5]. 

An approximation of the particle's position X t is obtained in [S] to be X t ~ 
fit + \ftaZ as f ^ oo, where [i = k\v/{k2 + fci),cr = 2k2k\v 2 /{k2 + fci) 3 and Z is 
a standard Normal variable. That approximation is valid only for large fixed values 
of t, while we next extend this result to give an approximation for the whole time 
trajectory of the particle's path. This is accomplished by the following functional 
central limit theorem for the position of the particle undergoing stochastic transport. 

Proposition 3.4. Let X be the position of a particle following the piecewise 
linear Markov process from Proposition \3.3\ started at Xq — on transport, then 



in distribution on the space of continuous functions. 

Proof. For these results we use the notion of stochastic averaging [23] [24] [26] . 
Note that the indicator process £ for being on- or off- transport is independent of the 
position X of the particle. Hence, the position of the particle X is a linear random 
evolution process, Ch 12 of [TT], [TH], driven by the independent indicator process £. 
The generator of (X, £) is the closure of the operator 




and, if B denotes a standard Brownian motion 




Af{x, - o-{0d x fix, + A(£) [f{x, s{0) ~ f(x, 0] 
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for all / g V{A) = C i,u (R + x {0,1}) (the space of all continuously differcntiable 
functions in x continuous in £ and vanishing at infinity), where 

<t(0= v£, A(0 = fc 2 e + fci(l-0, and = 1 - £ 

(when £ = 1: a = v, X = k 2 , s = and when £ = 0: a = 0, A = fci, s = 1). In other 
words, if (Xo,£o) = (0, 1) we have 

x t = j\t.ds, & = i(i + (-1)^/0* 

where Y is a rate 1 Poisson process, and Y(f* X(^ s )ds) is a counting process of the 
number of switches of £ until time t. 

Rescaling time and position of the process by 1/n, we get that (^-, £ n -) satisfies 



X 



nt 

n 



J*vUds, U = ^(1 + (-l)^"^ 



and its generator is 

A n /(x, - v& x f{x, + n(fc 2 £ + (h (1 - 0)) [f{x, s(0)- f^, 0] 

for all / e £>(A") = Cq'°(M + x {0,1}). Note that £„. switches at rate propor- 
tional to n, forming an ergodic Markov chain with stationary distribution 7r(l) = 
fci/(fc 2 + fei),7r(0) = fc 2 /(fc 2 + fci), and / a(£)n(£) = vf £tt(£) = ufci/(fc 2 + fci). Hence, 
the strong ergodic theorem implies that 

x n i r . , vki 

— = - / v£ s ds — > - — — a.s. 
n n J n-s-oo fc 2 + k\ 

We can extend this to a functional statement on any finite time interval [0,f]. Fix 
t > and take any A > 0, then there exists jt-a < oo a.s. such that 

\X n vki | A 

< — , for Vn > nA 



fc 2 + fci t ' 



Now, let M = sup n< „ A |X„ — fci/(fc 2 + fci)tm|, which is finite a.s. since riA < oo a.s. 
Let n > M/A. Then, for any < s < t we have that either ns > n\ in which case 

vki i \(X ns vki \ I A 

s < — s < A 



1 n fc 2 + ki 1 \ ns fc 2 + fci / 1 i 

or, ns < riA in which case 

I X ns vki I 1 1 „ wfci , M 

1 — n~ s = - \ x ns - — n-nsl < — < A 

1 n fc 2 + fci 1 n 1 fc 2 + fci 1 n 

implying that we have sup 0<s<t |^;^ns — vki/(k 2 + fci)s| < A whenever n > M/A. 

Once we rescale the position for the particle by l/\/n and time by 1/n, £ still 
changes at a much faster rate than the position of the particle X. The generator of 
the rescaled centered process (X n ,£ n ) defined as 
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is the closure of the operator 



)d x f(x,0+nX(0(f(x,s(0) - f(x,0) 



on / € T>(A n ) = Cq' (M x {0,1}). We will use the stochastic averaging theorem 
(Theorem 2.1 of [26]) to show that the paths of the centered rescaled process converge 
in distribution to paths of a Brownian motion with a diffusion coefficient equal to 
2k 2 k 1 /(k 2 + k 3 ) 3 v 2 . 

Let h(£) be the function 



HO = v 



k 2 ki 



1 k 2 ki 
= v- 



1 



(fc 2 + fcl) 3 A(£) (*2 + *l) 3 + (i-0 



Then h(l) = vk 1 /(k 2 + fci) 2 , h(0) = -vk 2 /{k 2 + hf imply that h(s(l)) - h(l) = 
-v/(k 2 + kt), h(s(0))-h(0) = v/(k 2 + fci), which in turn imply that A(l)(/i(s(l)) - 
h(l)) = -vk 2 /(k 2 + ki), A(0)(/i(s(0)) - h(0)) = vki/(fo + fci), so that for £ e {0, 1} 



A(0(M«(0)-M0) = -K-« 



k 2 + ki ' 



Now, for any / G C 2 (R) define a sequence of functions f n G Cq'°(R x {0, 1}) by 



f n (x,Z) = f(x) + ^h(Od x f(x) 



Then /" — > / as n — > oo and 



k 2 + fci 

^A(0(ft(*(0)-fc(0)Qr/(a: 
fci 



(v£ — v 



k 2 + fci 

where A is defined on V(A) = C^(R) by 



)/i(0^/(a=) = m*) 



k 2 k\v 2 
(k 2 + fci) : 



r^/(a 



Define a sequence of processes 



1 



— h(tf)d x f(x?) = f n {x?,%) - /(*») 



Then our earlier calculation implies that for any / G T>{A) 



Af(X2)ds + e{' n = r(X?,$) 



A n f(X?,C)ds 



is a sequence of martingales. Since / G C 2 (R), £™ G {0, 1} it is clear that 
rt 



sup E 



< oo and E 



sup |e 

s<t 



f,n\ 
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In order to apply Theorem 2.1 of [55] on stochastic averaging it is only left to show 
the process X n satisfies the compact containment condition, that is for any t > and 
A > there exists a compact set K C M such that 

inf P{X™ €KVs<t}>l-A 

n 

This follows from the fact that X" + h(^)/^/n is a sequence of martingales (let 
f(x) = x) with mean 



E 



and second moment (let f(x) = x 2 ) 



E 



XI 



k 2 k\ v 2 
y/n (k 2 + kiYy/n 

v 2 2t+ nmA 



So, by Doob's inequality 



sup 

8<t 



x: 



X? 



(k-2 + ^f 

hm\ 2 



) c 2 ( 



k 2 ki 



v l 2t- 



Noting that h m i n := v min(fc 2 , ki)/(k 2 + k\) 2 < h < h max := v max(fc 2 , ki)/(k 2 +/ci) 2 , 
and choosing C (given on t and A) so that the right hand side of the inequality with 
Ti = 1 is less than A, shows that with K = [— C — h max , C — /i m in] the compact 
containment condition holds for (I")„>i. 

Now Theorem 2.1 of [26] implies that X n =5> W in distribution on the Skorkhod 
space of continuous functions, where W is a process with generator A, and conse- 
quently has the same distribution as y '2k 2 kiv 2 / '(k 2 + ki) 3 B where B is a standard 
Brownian motion. □ 

3.3. Connection to Partial Differential Equations Models. In order to 
demonstrate the connection between our model and the PDEs seen in [36] | 35 j [IB] . 
consider the process (X, £) of the particle following the piecewise linear Markov pro- 
cess, and for any x > 0, t > let 

q(x, t) = P{X t e dx, £ t = l}/dx, p(x, t) = P{X t e dx, ^ = 0}/dx 

denote the probability densities of the particle's location x on and off transport, 
respectively, over time. Kolmogorov forward equations for (X, £) imply that q and p 
satisfy the system of PDEs 



d t q(x, t) + vd x q(x, t) = -k 2 q(x, t) + k\p{x, t) 
d t p{x,t) = k 2 q(x,t) - kip(x,t). 



(3.1) 
(3.2) 



When k 2 = = k\, the limiting PDE is simple linear transport: (dt+vd x )q(x, t) = 
0. The initial condition q(x, 0) = Sq(x) corresponds to the density of a single particle 
at the origin at t = 0. The time evolution via simple linear transport is translation of 
the delta function, while the time evolution via the equations p.l[) and (|3.2[) will have 
a spreading profile. This is clear from the macroscopic limits of (X s , £ 5 ) as 5 — > 0. 
When k 2 — ki — the particle never switches off from traveling on transport at speed 
v and is deterministic, as seen in Proposition [3H] When k 2 ,ki > the particle follows 
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a truely stochastic process (X, £) with a non-zero variance, as seen in Proposition 13.31 
and Proposition ^. 41 

In the experiments described in the introduction one sees "approximate" traveling 
waves of radioactivity in the axons in the sense that there is a slowly spreading wave 
front moving at constant velocity away from the cell body. Equations (|3.1[) and (]3.2[) 
are linear and do not have solutions that are bounded traveling waves. It was shown 
by a perturbation theory argument in [36] [35] that as e — > the solutions of 

e(d t + vd x )q s (x, t) = -k 2 q £ (x, t) + kip e (x, t), (3.3) 
ed t p £ (x, t) = k 2 q £ (x, t) - k lP £ (x, t). (3.4) 

subject to q £ (0 7 t) — qo, are to leading order 

q £ (x, t) = cxH( X ' ,t), p £ (x, t) = c 2 H( X ' ff ,t), 

where H satisfies the heat equation 

d s H(y, s) = —d yy H{y, s), H(y, 0) = x ( _oo,o) (3-5) 



k 2 v 2 2k 2 k\v 2 k\ k 2 

^ k 2 + ki a {k-z + kiY Cl k 2 + k\ ' ° 2 k 2 + k x ' 

This asymptotic form is valid for small e, that is for large k 2 and k\. However, if 
we set q(x,t) = q e (§, f) and p(x,t) = p £ (j, f), then q and p satisfy ([3~3]) and (f3T4|l . 
so the solutions of (13. ip and (|3.2I) behave like approximate traveling waves for large 
t and large x whether or not k 2 and k± are large. These results have been proven 
rigorously by Friedman and coworkers |14] |15) [16) [17] . 

To see that our Proposition ^. 4l provides another rigorous proof of these properties, 
albeit using stochastic methods, consider the process (eX./ E ,^./ E ) (that is (jX„.,^„.) 
with n = 1/e), and let 

q £ (x,t)=P{eX t/£ g dx,&/e = l}/dx, P £ (x,t) =¥{eX t/e G dx, £ t/e = 0}/dx 

be the probability densities for this process. The generator of this process is A n 
(n = 1/e), so the Kolmogorov forward equations imply that q e and p £ satisfy the 
system of PDEs p.3[) and (|3.4p . Our result from Proposition ^. 41 states that 

P{eX t/E G dx}/dx a H(^J^,t) 

for small e > 0, where if satisfies ([33]) . Hence, q £ +p £ w £f(^/r>*)> andP{£ f/e = 1} « 
ki/(k 2 + fci), P{Ct/e = 0} ~ k 2 /(k 2 + ki), gives the result that g e (cc, i) and p e (x, i) are 
well approximated by ^-H(^,t) and ^M_H(^,t). 



4. Dynamics from the spatial system perspective. 

4.1. The spatial system in equilibrium. We are now ready to characterize 
the steady state dynamics induced by continually adding particles from the nucleus 
and removing them when they reach the distal end of the cell. 

Proposition 4.1. Let ({(Qi(t),Pi(t)),i = l,...,N}) t >o be the number of parti- 
cles in the axonal transport system with compartments of size 5, on and off transport. 
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respectively. Suppose the rate of production of particles from the source is rq a — vp . 
Then this Markov chain has the product-form stationary distribution 



{Qi, Pi) ~ Pois(q ) <g> Pois 



fc2(?o 



where all {(Qi, Pi),i = 1, . . . , AT} are mutually independent. 

Proof. Since the production rate is rq , the generator of the process (Q\,P\) is 

A qo f(q,p) = [f(q + l,p) - f(q,p)]rq + [f(q - l,p) - f(q,p)]rq 

+ [f(q - l,p+ 1) - f(q,p)]k 2 q + [f(q + l,p- 1) - f(q,p)]k lP 



If we use f(q,p) = Qi(t) and f(q,p) = Pi(t), and take expectations, we get a system 
of ODEs governing the change in E[Qi], E[Pi] over time 



<IE Q 1 ] ( l]) = rq Q - rE[Q a (i)] - AfcE^t)] + fciE[Pi (t)] 



dt 

dEjP^t) 
dt 



= h 2 E[Q 1 {t)]-k 1 E[P 1 {t)] 



indicating that in equilibrium in the first section the mean numbers of on-transport 
particles and off-transport particles are E[Qi] = go, and E[Pi] — E[Qi]|^- = go^Ai, 
respectively. Let ir qo (g, p) — tt\ q (q) <g> n\ p (p) be a product of two independent Poisson 
distributions with rates Ag = qo, an d Ap = gofe/fci respectively. To show that 
ir qo (q,p) is a stationary distribution for the process (Qi, Pi) we need to check that 



00 00 



q=0 p=0 



for any choice of function / e V(A qo ). 



00 00 \ij ,p 



00 00 



~ (Aq+Ap) EE 44( [/(<z + 1,p) ~ + tffa - ^ - /fo^to 

+ [f(q - l,p+ 1) - f(q,p)]k 2 q + [f(q + l,p- 1) - /( 9) p)]fcip) 
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In the two sums the factor multiplying /(g,p) for any (g,p) € N x N comes only from 
terms involving {q — 1, q, q + 1} and {p — l,p,p + 1} and equals e - ( A<5+Ap ' times 

Aq 1 A p A Q A P A^ +1 A p A Q \ p 

-rq t-T t 1o + 7 — — 777 + : ) 1 ~ 1 r 9 



(g-l)!p! 9! p! (g + l)!p! v g! p! 

-k 2 {q + l) f-T k ^1 



(g + l)!(p-l)! ^ ' g! p! 

rfci(p+l)- -T-T-kip 



(g-l)!(p + l)! — ' g! „! 
g! p! V Aq g + 1 g + 1 Ap 

- fc2<?+ . 9 Ap - fci(p + i) - hp 

Agp + 1 



A^A|/^ 
g! p! Vg 



rgo - rg + go?" - t-?) =0 



since Aq — go and Ap/Ag = kz/ki. 

Thus, in equilibrium the input rate for (Q2,-£2), which is rQ\, has a Poisson 
distribution with mean rgo , and is independent of P\ . Let 7r go (gi ) ® ir\ Q (q 2 ) ® 7r,\ p (P2 ) 
be a product of three Poisson distributions with rates go, Aq = go, and Ap = gofe/fci 
respectively. To show that this is a stationary distribution for the process (Qi, Q2, P2) 
we need to check that 



00 00 00 



Y Y Y Ai/te^KoteKite^) = 

91=0 92=0 P2=0 

for any choice of function / e 0(^4 gi ). For each hxed value gi, according to our 
previous calculation the inner two sums give 0, so the whole sum is 0. 

Thus, in equilibrium, {Q2,P 2 ) have the distribution it\ Q (q) ® t^x p (p) with Aq = 
90;Ap = qok 2 /ki as well, and are independent from (Qi,Pi). It follows by induc- 
tion that the stationary distributions for {(Qi,Pi)} are independent and identically 
distributed as ir\ Q (q) <g) ir\ p (p), with Aq = go, Ap = qak 2 /ki. This is also an exam- 
ple of a clustering process satisfying the detailed balance conditions with linear rates 
discussed in Sec. 8.2 of [25]. □ 

We point out that the mean number of particles both on and off transport is 
(1 + |j-)go = (1 + fa)po$: scaling with the size of a compartment. To obtain the mean 
number of particles per unit length we add particles in ss 1/5 compartments, and the 
mean number of particles per unit length is (1 + jr)po independent of the choice of 
compartment size. 

One immediate consequence is the analogous result for the number of particles 
in the stochastic non-compartmental model at any location along the axon. Namely, 
suppose the particles move according to the piecewise-linear Markov process (X, £) 
with a Poisson rate po influx of new particles at location 0. Then, at any location 
< x < L along the axon, the numbers of particles (Q( x ,x+dx)> P(x,x+dx)) on an d 
off transport, respectively, have the stationary distribution Pois(podx) (E)Pois(^ podx) 
where for any x\,...Xk G (0, L), {Q Xi }i<i<k and {P Xi }i<i<k are mutually indepen- 
dent. We note that this result would not have been obvious to notice without going 
through the compartmental model first, yet its consequences for prediction and anal- 
ysis of the long term stochastic behavior of the system are quite powerful. 
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4.2. Homogeneity of the axons at equilibrium. Recall that 6 = lOnm, 
roughly the step size of motor proteins, and that axons can be up to one meter in 
length. Thus we are interested in phenomena on all the length scales lO^J, where 
v = 0, 1, 2, ... 8. Let A = 10"8; we want to determine how similar different segments 
of the axon of size A are. Let Qa and Pa denote the numbers of on-track and off-track 
particles in a segment of length A. 

In equilibrium, Qa and Pa are both sums of 10" independent Poisson random 
variables with parameters Xq — qo and Xp = j^-qo, respectively. Therefore the distri- 
butions of Qa and Pa are Poisson with parameters 10 v Xq and 10 u Xp, respectively. 
The mean and the variance of the number of particles in the segment of length A is 
IO^Aq + Xp). To see how homogeneous different slices of length A are, we consider 
the coefficient of variation, ca , which is the standard deviation divided by the mean. 

1 1 

ca = — , = — — 

^(A Q + A P )10- y/(l + fc 2 /fci)<Zol0 1 ' 

As indicated in Section[2] qo is m the range 0.25 to 2.5 in different axons. For illus- 
trative purposes here, we will assume qo = 1- Since fe/^i = 3, we see that the scale- 
dependent coefficient of variation ca = l/(2v / 10 8 A). Therefore at the ten nanometer 
scale the coefficient of variation is simply 1/2. At the micron scale ca = 1/20 and at 
the millimeter scale ca = 0.5 x 10 -5 / 2 . The cutoff between "high variance" and "low 
variance" distributions is usually considered to be when the coefficient of variation is 
near 1, so by this standard the axon is extremely homogeneous in its length at large 
scales. 

4.3. Balance between efficient transport and targeted delivery. The pre- 
ceding characterization of the transport apparatus enables us to address questions 
concerning whole cell function. One core issue is that intracellular transport must 
simultaneously accommodate two functional demands: some material, such as the 
enzymes used to construct neurotransmitters, must be transported from the soma to 
the axon terminal in a timely manner; whereas other cargo, such as sodium channels, 
need to be delivered to unspecified locations as needs arise throughout the length of 
the cell. The tradeoff between these two goals is clear. If a typical vesicle spends the 
vast majority of its time in transport mode, the mean velocity will be close to the 
mean on-transport velocity, but any needs that arise in the central part of the cell 
will be neglected. On the other hand, if a typical vesicle spends too much time off 
transport, presumably available for use if needed locally, then it will take substantially 
longer to traverse the entire cell. 

Recently Bressloff and Newby [311 13"2"] modeled particles that are created near 
the nucleus that then undergo intermittent search (being in search mode while off 
transport and not in search mode when on transport) for a target hidden somewhere 
along the axon. Their model is a non-compartmental individual particle model with 
the additional feature that vesicles can move backwards as well as forwards. They 
compute the probability that the particle is successful and conditioning on success, 
the mean first hitting time. With our system-wide model we can accommodate the 
observation the if a given vesicle misses the target, another vesicle with similar cargo 
will pass by before too long. We will assess this hitting time under two assumptions 
about the density of relevant material. Our standing assumption qo ~ 1 is appropriate 
for types of cargo that are found densely throughout the cell. In this setting, the wait 
time is essentially just the time it takes for one of several nearby vesicles to unbind 
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from transport in the target region. A more interesting case is a setting where the 
needed cargo is sparsely distributed, say qg ~ 10~ 3 . In this setting, if the first cargo 
to reach the target region fails to unbind, there will be significant time before the 
next arrival. As we will see, we can still assess the trade-off intrinsic between risking 
a target miss and diminishing the time of the next arrival. 

To make the discussion concrete we define a target region R n — {i* + l, . . . ,i*+n} 
where u G 1, . . . , N — n. At time zero, we take the system {Qi(Q), -Pj(O)} to be drawn 
from the stationary distribution described by Proposition ^. 1 1 conditioned on the event 
that for all i € R n , Qi(0) = Pi(0) = 0. We introduce the hitting time H n := inf{£ > 
: YlieR ^»(*) > 0}, which marks the first time a particle is off transport while in 
R n . Since computing the mean of H n is analytically intractable, we introduce another 
hitting time H' nl stochastically dominating H n , that nevertheless reflects the essential 
tradeoff between maximizing mean velocity and making detachment from transport 
likely in the target region. 

Let /* :={ie {1, . . . , z*} | Qi(0) +Pj(0) > 0} be the set of all non-empty sections 
of the cell at time 0. Among the particle in these sections some will be "successful" 
in that they will detach from transport in the target region, while others will be 
unsuccessful. Labeling each particle in these terms, we decompose Iq into locations 
with successful particles I* and unsuccessful particles I". We are interested in the time 
H' n at which the rightmost successful particle, that is, a particle starting from position 
i m = max{/|}, detaches while in R n . If S = then we define i m = max{0} := 0. 
For the position of this particle we use the notation (A t 5 ) t > from Section [31 where 
Aq = Si m , X? € {5i m ,...,5N}, together with the indicator (£f) t >o,€t G {M} of 
whether the particle is on or off transport, respectively. Let 

H' n := M{t > : (A t 5 ,^) = (*(*. + 1), 1)}. 

Note that H' n = inf{t > : G 5R n }, since particles cannot skip sections, and 
always enter a section on transport. If we were to analogously define a sequence of 
times {H^} ieI , where for each i e U\ = inf{t > : X$ e 5R n } with Aq = Si 
and £q € {0,1}, then the exact first hitting time of the target region will satisfy 
H n = mm{#; : i € I*}. Hence, clearly H n < H%" = H' n . 

Since the exact distribution of H % n is complicated, E[i?„] is analytically intractable, 
and instead we focus on finding a simple expression for E[Z/^J. Our point is that, in 
the sparse material limit (qo small), the rightmost particle becomes increasingly likely 
to be the first successful particle to detach in the target region, and H n approaches 

The computation of E[iJ'J requires computing the time it takes a particle to travel 

a certain axonal distance, given by the following. 

Lemma 4.2. Let the initial position of the particle be (Aq,£q) = (0,1) and let 
G {1,...,L} be given where L is the total length of the axon. Then, the time 

T Lm = w£{t > : Xf = L*}, satisfies 

Proof. Since (A t 5 ,£j?) t >o is a Markov process with generator 

A 5 f(x,0 = rt[f(x + 5,Z)-f(x,0] 

+ k 2 t [/(*, t - 1) - f(x, 0] + fci(l - [f(x, £ + 1) - f(x, 0] 
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it follows that Mj 



rS J* gda and M ( 2 := # + /* M 5 s ds - J* fe x (l - g)da are 
both martingales. Since both M 1 and M 2 have bounded increments and E[TlJ < oo, 
the optional stopping theorem implies that 



= E[M 1 ] =E 






r f TL * 




r r Th - i 




= L* - vE 


/ £ds 


E 


/ Z s .ds 








-Jo 




-Jo J 



and 



i = e[m 2 ; 



E 



1 - fciEpYj + (fc 2 + fei)E 



which implies 



E[T L 



ki 



E 



and our claim follows. □ 

We also note that the same computation holds for the mean time a particle in 
the stochastic non-compartmental model {Xt,£t)t>o takes to reach a distance L*, as 
the two martingales used in the proof depend only on v and not on S. 

We next compute the hitting time H' n of a particle started at location i m G /* . 

Lemma 4.3. Let the system {(Qi(0), Pj(0)), i = 1,...,N} have the stationary 
distribution given by Proposition ^- 1\ conditional on Qi(0) — -Pj(O) = 0, for all i G R n - 
Then 



(i 



-An 



1 



rqoPn k 1 (k 2 + k 1 ) 



kip-n 



k 2 



1 



nk 2 
k 2 + r 



where p n 



1 



■k 2 +r 



)™ and X n 



— — Qo l *Pn- 



Proof. The proof of Proposition 14.11 shows that conditioning on the values of 
Pi(0),i G BP does not affect the law of {{Qi(t),Pi(t)),i < i*}, hence for any t > 0, 
they form two mutually independent sequences of Pois(qo) a- n d Pois(-^|Y e ') random 
variables. 

Let S denote the number of successful particles between sites and i* at time 
zero. The total of number of particles at time at these sites that are either on 
or off transport is distributed as a Pois( fc2 ^ fcl 9o»*) variable. We next compute the 
probability p n that any particle once it reached the target region is "successful" in 
detaching there. This probability can be written p n := Yn=i P(^) wnere p{i) 1S the 



-1 k 2 



is the chance 



< r y 

\k 2 +r' k 2 +r 

= E^Ki) = i-(^F) T 



probability it first detaches at location i» +i. Since p(i) = 
a particle gets off transport in the i-th compartment, p, 

The probability that a particle is successful does not depend on which location it 
was at time 0, and whether it was on or off transport at that time. Hence, S is is 
distributed as a Pois( fc2 ^ fcl qoi*p n ) variable, and conditioned on the value S, the set 
of locations of the particles at time is distributed as a set of S draws from the 
uniform distribution on {1, . . . Since i m is the maximum of S samples from a 

Uniformjl, . . . , i»} distribution we have E[i m |S] = i* — J2 % x=i xS ~ s+T- 

We now decompose H' n = H' e n + H' o n , where we let H' e n is the time it takes a 
successful particle initially at location i m to enter the region i? n , and H Dtn is the time 
it takes any successful particle after it enters the region to get off transport. Thus, 



[H'n\ 



E l H 'o,n\ 
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If a successful particle starts at location Xq 



£o=l, then the time it takes to enter the region R n is by Lemma 



G {0, ...,&*} on transport 



E[ J ff;j(x 5 ,^) = (* m ,i)] 



(i* - i m )8(k 2 + fei) 



If it starts at location Xq — i rn e {1, . . . , i*} off transport £q = 1 then it takes an 
additional Exponential time with mean 1/fci for it to get back on transport at the 
same location, so E[H' e J(X*, = (i m ,0)] = E[H' ein \(X(>, $) = (i m ,l)] + and 
since we assume new particles always enter the system on transport 



E[<nl*'« 



(i* - i m )S{k 2 + fci) 



7— li m >0- 

fci 



Since E[i m \S] = Ugfj, F{i m > 0} = P{5 > 0}, we get that 



EK„] = 

with 5 ~ Pois(A„), A„ = 
we get 



S + l 



j»(fc 2 + fci) 
rfci 



fci(fc 2 + fci) 



fc2+fcl 



^q i*p„. Since E[5] 



(1 



PIS' > 0} 



B {5 > 0} = 1 - e 



-An 



A* 



rgoPn fci(fc 2 + fci) 



To calculate E [i?^ n J note that if a particle first gets off transport in the £-th 
compartment then the time of its travel until this point is a sum of i independent 
and identically distributed exponential random variables with parameter fc 2 + r. The 
probability a successful particle first gets off transport in the z-th compartment is 
Hence, the time a successful particle takes to get off transport once it enters the 
region R has the mean 



i=l 



Pn 



nki 



and our claim follows. □ 

To complete our analysis, we wish to characterize this result in terms of length 
along the axon and independent of the stepping size 5. To this end, we fix a length I 
and for a given compartment size £, we let n = [^/<5], and for the start of the region 
we let z* = \£*/5] . As such, as 5 — > the limiting region becomes Ri — (£*,£* + £). 
Then, under the assumption that qq/5 — > po we have 



and E[H' n \ -> E[H'\ where 



-£k 2 /r 



A-n — > Xf 



k 2 + fci 
fci 



Po£*Pe 



E[H'] = (1 - 



1 



vpope fci(fc 2 + fci) 



e-^(l + ^) 



(4.1) 



It remains to interpret this result with respect to the parameter choices we have 
made. First, we must set a value for the size I of the target region. For this purpose 
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we note that the typical size of a Node of Ranvier - a gap in the myelin sheath of a 
myelinated axon where sodium channels are concentrated in the cell membrane - is 
approximately one micron^ 

As seen in the preceding proof there are three contributions to E[iJ'], and we next 
analyze them in terms of their dominance for the overall value. We begin with the 
contribution from the last term, which is the time it takes for a successful particle to 
detach once it has reached the target region. By our choice of t, the recurring ratio 
tjv is one. Along with the earlier assumption that k 2 = 1, the entire term simplifies 
to (e — 2)/(e — 1) = 0.4 s. The multiplicative factor 1 — e~ Xe preceding the first term of 
(14. ip results from a boundary effect: when is very close to zero, it is very unlikely 
there are any particles already in the system between the soma and the target region. 
When the target region is in the middle of the axon, this contribution from particles 
in that section of the axon at time is significant. 

When qo ~ 1 as assumed earlier, then po = qo/6 ~ 10 8 . Then if t„ > 10~ 8 , which 
corresponds to the target region being just one motor step down the length of the 
axon, we have po^* ~ 1 an d > 2.3 so 1 — e _Af > 0.9. Looking at the first term 
inside the parenthesis, we note that vp ~ 10 2 , while pg — 1 — ex-p(—£k 2 /v) « 0.6 is 
the probability that any given vesicle will be successful in detaching from transport in 
the target region, so l/(vpopi) ~ 10~ 2 . Meanwhile, the term ^2/(^1(^2 + which 
is the expected time to bind to transport if a successful particle happens to be off 
transport at time 0, is 1/4 for the assumed values of k 2 and k±. Therefore under the 
go ~ 1 assumption, both this and the final term contribute significantly to the hitting 
time. 

However, in the sparse material regime, say go ~ 10~ 3 , we have then po ~ 10 5 , 
and the factor 1 — e~ A * > 0.9 when I* > 10~ 5 . That is, if the target region is at 
least 1/100, or at least 1000 segments, down the length of a 10~ 3 ?n axon. However, 
now the rate limiting factor is the wait time for the first successful vesicle to arrive 
in the target region, which is captured by the term 1/vpop. The product vpo « 0.1 
measures the average rate at which new particles should arrive and together with the 
probability of success pe ~ 0.6, we now have 1 / vpope « 16s. So, in the sparse material 
regime the first summand in E[i?'] giving the mean time for arrival of the particles to 
the target region dominates. 

It is interesting to consider what happens to the arrival rate term under per- 
turbations of the various parameters. In particular, we draw the reader's attention 
to changes in v. When viewing intermittent search in terms of a single particle, the 
probability of finding the target is strictly decreasing in v. Higher velocity seems to be 
the enemy of finding the target. Indeed, from a system point of view, this implies that 
the system will require more trials before a successful particle arrives in the target 
region. What the equation (|4.1j) gives us is the ability assess how much more quickly 
the trials will happen. In fact, the function v(l — exp(—£k 2 /v)) is increasing in v. 
Therefore the entire expected wait time E[ff'] is actually decreasing in v. That is to 
say, while the particles are less likely to succeed, they will be arriving more rapidly 
enough to counterbalance the lost time. We believe that this kind of quantitative 
analysis will prove fruitful in future study when coupled with more details of the 



1 Wc do not wish to claim this is a complete model for deposition of sodium channels in Nodes of 
Ranvier, as the particulars of the biology - which may include factors like local signals that encourage 
motors to detach from microtubules near the Nodes - are not fully understood. We merely wish to 
use the size of the Nodes to fix our intuition about the size of a target to other important length 
scales such as that of microtubules which are also a micron in length. 



21 



biology of deposition of materials in the cell membrane. 

4.4. Approaching Equilibrium. We have seen above that the axon is very 
homogeneous at stochastic equilibrium on a space scale down to micrometers. One 
of the beautiful properties of transport with reversible binding is that if it is locally 
out of equilibrium, the on-off dynamics can return the system to equilibrium on a 
much faster time scale than waiting for new material to arrive from the nucleus. 
Furthermore, as is discussed in [32) , one of the key goals of biophysics investigation is 
the discovery of behaviors for which biochemical regulation is not necessary. This is 
of fundamental importance to the biological function of the system because it means 
that the axon will automatically "repair" itself without central control of the repair 
process. 

How good is this mechanism? If a segment of the axon is far away from equilib- 
rium, how long does it take to get back close to equilibrium? To investigate this ques- 
tion, we imagine that the axon is at stochastic equilibrium except for some segment 
R, where the total number of particles on and off transport, Qj(0)+Pj(0) = 0, Vi G R, 
is zero at time 0. Let a and e be given small numbers and suppose that is the 
mean vector for (Qi,Pi) at stochasttic equilibrium. We want to compute (an upper 
bound for) the time t* so that 



that is, the probability that (Qj,Pj),Vi G R is significantly different from is very 
small. In the applications below we will choose a = 0.1 and e = 0.05, and we will 
see that a 10 micron segment can recover in about 10 seconds, while a 1 millimeter 
segment will take about 1000 seconds or 15 minutes to recover. We study this question 
first for a single location R = {i}, and then use the estimates derived to scale the 
results to segments of any length. 

PROPOSITION 4.4. Let (Qi(0), P»(0)) = (0,0) and let the constants a > and 
e G (0, 1) satisfy the relationship a 2 e|A 00 | > 1 where A^ = qo(l, j^) is the equilibrium 
vector of (Qi,Pi)- Then there exists t* > such that Vi > f*, 



We note that given a particular choice of a and e, the condition a 2 e|Aoo| > 1 
guarantees that there are "enough" particles. 

Proof. Note that for any given j3 G (0, 1), we may choose > such that for all 
t > t*, the vector of means X(t) := E[(Qj(i), Pi(t))] satisfies |A(t) — Aoo| < a/3 1 Aoo |- 
Then 



F{\(Qi{t),Pi(t)) - A^l > alAool} < F{\{Qi(t),Pi(t)) - X(t)\ + \\(t) - A^ > a^} 

< P{|(Q i (t),P i (t)) - \(t)\ > a(l - P)\Xoo\} 



Applying Chebyshev's Inequality, and observing that the variance of a Poisson random 
variable is equal to its mean, we conclude that 



P{|(Qi(t*),Pi(t*)) - Xoo\ > olAool.Vi G R} < e, 



P{|(Qi(t),Pi(*))- Aool > a\X x \}<e. 




a : 




P{|(<3i(t),P(*))-A oo |>a|A 0O |} < 



Var[|(Q t (t),P(t))|] 
a 2 (l-/3)2|Aoo| 2 



a 



2 (1 -/3) 2 |Aoo| 2 



|A(t)| 
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Since the initial condition for both P$ and Qi are less than their respective equilibrium 
values, each are monotonically increasing functions and the above reduces to 

1 



P{|(Q l (t),P I (t))-A 00 |>a|A 00 |}< 



a 2(i-/3)2| Aoo | 

for all t > t*. In order to satisfy the requirement that the right hand side must be less 
than e, we solve for j3 and find j3 — 1 — (ay^F|A^|) _1 provided that a-^/elAool > 1. 

It remains to study the convergence of the mean and the appropriate choice of 
t*. The dynamics of the mean vector A(t) are given by the ODE 

j t \(t) = -AiA(i) + gorei (4.2) 



k 2 + r —ki 
-k 2 ki 

The solution to (|4~2|) is A(i) = X x + e~ Alt (\(0) - Xoo) where Aoo = q rA^ 1 ei = 



where ei is the unit vector (1, 0) and A\ 

n 

qo(l, I 2 -)- This yields the estimate 



where a is the smaller of the eigenvalues of A\. Noting that a > 0, t» may be chosen 
so that e~ at * < a/3, i.e. t„ = a -1 ln(l/(a/3)) is an upper bound for the time to be 
close to equilibrium with high probability. □ 

In order to calculate return to equilibrium at various scales, we now suppose 
that the whole axon is in statistical equilibrium except for a segment R of length 
A = S10 1 ' in which we will assume that there are no particles either on or off the 
track. Proposition 14.41 covered the case v = 0. We are interested in v = 1, . . . , 8. We 
imagine that the axon is broken up into 10 8_!/ segments of length A. In this rescaled 
system, the unbinding and binding rates per particle, k 2 and k\ remain the same, as 
well as the mean on-transport velocity v. In order to retain this mean velocity, the 
rate of lateral stepping must be decreased to r — rl0 _ly . 

The ODE for the mean vector of the rescaled system is given by 

^A(t) = -IiA(i) + g rei (4.3) 
at 



with Ai 



fc 2 + r —k\ 
~k 2 k\ 

We note that the last term in (|4.3[) contains an r rather than an f. This is because 
the input rate is unchanged while the exit rate is diminished. The resulting equilibrium 
value is therefore rescaled as well, A,^ = q^rA^ei = 9o^( fe2 / fcl ) = 9olO I/ ( fc2 J fcl ) ■ Both 
components of this vector are of order 10^, as expected. 

Using the parameters discussed in Section [21 k 2 ~ 1, k\ = = 10~ 6 m/s,r = 
10 2 s _1 . This implies f = 10 2_I/ s _1 . We choose the thresholds to be a = 0.1 and 
e = 0.05, so the constraint that a 2 ^^^! > 1 requires that |Aoo| > 2 x 10 3 . Since 
0.25 < qo < 2.5, this is indeed the case if v > 3, i.e. if the segment has length greater 
than 10 microns. 

It follows from Proposition 14.41 that the time to equilibrium, t* , is proportional 
to a~ 1 where a satisfies: 

a = - {k 2 + fei + f - ^(k 2 + fci + f ) 2 - Afar) 

_ 2kir 

k 2 + k 1 +f+ ^{k-i + kx +r) 2 -4fcif' 
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For the given parameter values (with jA^I > 5 x 10 3 in particular), the constant 



of proportionality In ^•\/e|A 00 |/(a'\/e|A 00 | — l)j is contained in the interval (2.3,2.5) 

and does not have a impact on how the relaxation time scales with v. In terms 
of analyzing a, we note that kif is small compared to &2- To leading order, a ~ 
(fcif)/(fc 2 + ki + f) ~ 10 2-I/ s -1 . It follows that t* ~ 10"- 2 s. Thus, for a 10 micron 
segment (f = 3) the time to recover is about 10 seconds and for a 1 millimeter seg- 
ment (y — 5) the time to recover is 1000 seconds or 15 minutes. The time to recover 
depends, of course, on the parameter e that represents what we mean by "close." We 
also note that one can compute various measures of time to recover using the PDE 
models discussed in Section 3. 

5. Discussion. In this paper, we created a spatial Markov chain model for study- 
ing various aspects of fast axonal transport. Previous models that use PDEs treat 
the velocity of transport as constant when particles are attached to the fast transport 
system. Since it is known that transport along the microtubules is itself stochastic, 
it is important to have a fully stochastic model. Our model allows us to unify and 
extend previous work. In Section 3.2 we show that from the particle perspective as 
the compartment size S tends to zero, our model converges in distribution on the 
Skorokhod space of cadlag functions to the piecewise-deterministic model analyzed 
by Brooks [5]. Namely, we show that the paths of particles in our model converge to 
those of particles in a stochastic non-compartmental model. The argument proceeds 
by an explicit computation of the finite-dimensional distributions and a tightness ar- 
gument. In Proposition 13.41 we give a rigorous probabilistic proof of why the paths of 
particles follow "approximate traveling waves" described by other authors [36l [14] . 
This proof is based on stochastic averaging arguments which show that a functional 
central limit theorem holds on the space of continuous functions for the paths of 
particles as the compartment size decreases. The diffusion of particles around their 
mean position can consequently be approximately described jointly for all time by a 
Brownian motion with the appropriate diffusion coefficient. 

In Section 2, we show how to use existing experimental data to indentify (ranges 
for) all the parameters of our model. In light of this, we can use the model to 
investigate several important biological questions. These are based on describing the 
spatial distribution of multiple particles in our model. In Section 4.1 we derive the 
stationary distribution for the number of particles in different compartments on and 
off transport along the axon. This gives an explicit description of the stochasticity of 
the system that is present even after a long time. In Section 4.2 we derive estimates 
for how homogeneous the axon is on different spatial scales. In Section 4.3 we study a 
question introduced by Bressloff, by providing a stochastic quantity which describes 
the balance the system needs to achieve between rapid transport that brings new 
material quickly and efficient local search that improves time of delivery to a target. 
Finally, in Section 4.4, we use the model to calculate the length of time that it would 
take for axonal segments of different lengths to recover to near stochastic equilibrium 
after they have been depleted of vesicles. 

In our stochastic compartmcntal model all event wait times are assumed to be 
exponential random variables, but this is certainly a simplification. As an example, 
the stepping process of kinesin is a well-studied though still not completely understood 
phenomenon. Much work has focused on assessing the dependence of the mean rate of 
translocation on both the load and the local concentration of ATP [U] [35] . Implicit 
in this analysis is the assumption of exponential wait times with state dependent rate 



24 



parameters. However, when fitting to data and matching dispersion information the 
authors in [T3] found it necessary to generalize the wait time distribution. This was 
followed by more detailed models for which it was shown that load carrying could 
in fact regularize the stepping times of kinesin motors |37) [9]. Generalizing waiting 
times would significantly affect our results. Since the particle position process is no 
longer Markov, we no longer have the direct connection to the previous results stated 
in Section 3.1.-3.3., nor can we use the stationary distribution employed in Section 
4.1 and used for addressing the biologicals questions in Sections 4.2.-4.4. In light of 
the known need for generalized wait times in the stepping process, it seems likely that 
detailed observation of the rebinding process will call for new mathematical models as 
well. Recall that when a vesicle unbinds from a microtubule it is unclear whether it 
typically rebinds to the same microtubule or if it explores the region significantly via 
diffusion before finding a different microtubule to bind to. In the latter case, a more 
appropriate model for rebinding time would be to solve some kind of first passage 
time problem and use that distribution for the rebinding wait. 

An important aspect of the biology of axonal transport is not included in the 
model presented here, namely the local deposition and eventual degradation of trans- 
ported materials. For example, sodium channels and sodium pumps are synthesized 
at the soma, transported down the axon and deposited in the axonal membrane, ether 
uniformly as in an unmyelinated axon or at the nodes of Ranvier in a myelinated axon. 
Channels and pumps are proteins with half- lives on the order of days to weeks. The 
present model can be extended to include a deposition compartment at each location, 
and, clearly, the processes of deposition and subsequent degradation will cause the 
mean number of particles both on and off transport to be monotone decreasing as 
one moves down the axon. How inhomogeneous this makes the axon will depend on 
the details of deposition and degradation rates. Our preliminary calculations indi- 
cate that long axons, such as the meter-long axons in human sciatic nerve, would be 
quite inhomogeneous. This is an important biological issue because it is controversial 
whether the machinery for protein synthesis (i.e. ribosomes) exist in axons [42] . We 
have also not included retrograde transport or the fact that some axons may have 
location-dependent unbinding rates [TO] . All of these issues will be the subject of 
future work. 
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